A Decade of Cave Drip Hydrographs Shows Spatial and Temporal Variability in Epikarst Storage and Recharge to Appalachian Karst Systems

Abstract

We conducted recession analyses on cave drip hydrographs from a 10-year record (2008–2018) of three drip monitoring stations within James Cave (Pulaski County, VA, USA) to examine differences in hydrologic characteristics of the epikarst and quantify the storage volume of the epikarst feeding the drips. We used two recession analysis methods (correlation and matching strip) to calculate recession coefficients for multiple hydrographs at each site. Results show subtle differences between the three drip sites, suggestive of spatial heterogeneity in permeability and storage in the overlying epikarst. Storage volume calculations show that during the recharge season, up to 95% of recharge through the epikarst to the cave occurs through rapid pathways (i.e., fractures), and 5% of recharge occurs through diffuse pathways (i.e., pores). However, during the recession period, recharge through rapid pathways in the epikarst decreases and occurs predominantly through diffuse flow. Combined, these results underscore the importance of documenting spatial and temporal characterization of drip rates and other recharge inputs into karst systems.

1. Introduction

Karst aquifers, characterized by soluble rocks such as limestone and dolostone, provide drinking water to 20–25% of the world’s population [1]. Effective management of karst aquifers requires accurate information on the recharge, or rate of replenishment, of groundwater to the aquifer [2]. Recharge to karst aquifers can occur from both diffuse flow (also called matrix, slow or base flow) and rapid flow (also called concentrated, conduit, or quick flow). Diffuse recharge is derived from precipitation or snowmelt that infiltrates into the pores of the soil and bedrock above a karst aquifer, which can then recharge the underlying aquifer through cave drips [3,4,5,6,7]. Rapid recharge can occur through larger fractures, sinkholes or sinking streams, providing highly focused recharge. The relative proportion of diffuse to rapid recharge reflects the distribution and connection of conduits, the permeability of soils and bedrock and the timing of water fluxes, which significantly affects the variability in both the quantity and quality of aquifer recharge in karst environments [1,8].

Estimating diffuse recharge to karst aquifers is particularly challenging due to the heterogeneity of the epikarst (the zone of soil and weathered bedrock above a karst aquifer), which can have high secondary porosity [1,8]. Complex distribution of fractures and conduits within the epikarst matrix results in spatially variable water flow, particularly with depth [4]. The shallow epikarst can transmit water first through high hydraulic conductivity fractures, but transmission of water downward may slow as the epikarst matrix creates a bottleneck where pressure from stored water transmits recharge to the karst aquifer via diffuse pathways [6,9,10].

A variety of approaches have been used to evaluate recharge mechanisms to karst systems including tracer and geochemical methods, rainfall-recharge methods, numerical models, statistical methods, and recession methods. Each method has advantages and disadvantages depending on the system under study. In karst aquifers with spring outlets, recharge can be investigated using measurements of discharge and water chemistry of springs (e.g., [8]). Tracer studies in karst aquifers are often used to delineate recharge source zones and the timing of recharge components [8,11,12,13]. Geochemical analysis of spring water has been used to infer properties of recharge in karst systems, including the function of the epikarst for storing recharge [14,15,16] and tracking the timing and components of recharge for karst aquifers [6,17]. Rainfall-recharge methods have been applied to quantify precipitation thresholds that initiate cave drips in Australia [18,19] and Virginia, USA [7], among other locations. Numerical models can be used to quantify the flux of water into and out of the conduit system of the karst aquifer [20,21] and also to differentiate recharge ages using chemical tracers [22]. Statistical analyses have also been used to examine recharge in karst. For example, time series models can provide insight into the relationship between precipitation and spring discharge [23].

Recession methods, which involve analysis of hydrographs, can provide information about the distribution of recharge (temporal and spatial) and storage properties of the surrounding aquifer [24,25,26,27,28]. Recession analysis has primarily been applied to surface water systems [29,30,31] but it also has been utilized for karst springs to examine temporal distribution of recharge, including both diffuse and rapid components [9,20,24,32,33,34,35]. More recently, recession analysis has been conducted on cave drips to derive storage parameters for a karst aquifer [5], showing that recession methods can be useful tools for extracting a variety of hydrologic characteristics in karst systems.

The complexity of karst aquifers requires both qualitative and quantitative analysis to examine how changes in annual precipitation, which are predicted to occur over the remainder of this century [36], will affect recharge to these aquifers. The existing body of research on cave drips generally utilizes hydrologic datasets on time scales which range from hours to months. Few studies have collected high frequency cave drip data over longer time periods due to equipment, data storage and access issues. Using data sets on the decadal scale or longer is beneficial in that they allow for long-term trends to be evaluated, particularly in sensitive settings such as karst terranes. Additionally, long-term data sets allow for comparison across spatial and temporal scales which serves to identify sources of change [17,20,23].

In this study, we used a 10-year dataset on cave drips and precipitation at James Cave, Virginia, to study spatial and temporal patterns of recharge to the cave and underlying aquifer. We conducted recession analyses on cave drip hydrographs to examine differences in hydrologic characteristics of the epikarst and quantify the storage volume of the epikarst.

2. Materials and Methods

2.1. Site Description and Instrumentation

James Cave (Pulaski County, Virginia, USA) lies within the Valley and Ridge physiographic province (Figure 1). James Cave is entirely within the Conococheague Formation (Supplementary Materials, Figure S1), a fine-grained bluish-gray limestone with interbeds of shale and varying amounts of dolomite [37]. The site is in a temperate climate zone with an approximate mean annual temperature of 12.8 °C [38]. On average, Pulaski County receives annual precipitation of 41.0 inches (1042 mm) [38]. Evapotranspiration, estimated using the Penman–Monteith equation, is highest in the summer months [39].

A portion of James Cave was instrumented in 2007–2008 with three cave drip monitoring sites (MS1, MS3, MS3; Figure 1). The three sites were assembled to measure drips from the cave ceiling (see Supplementary Materials, Figures S2 and S3 and [40] for pictures of the monitoring sites). The three drip sites were established in locations where the cave ceiling lies approximately 15 m beneath the land surface (MS1 and 2 are ~10 m below land surface; MS3 is ~7 m below land surface; [41]), which corresponds to the range of epikarst depth in many karst regions. The drip tarps were assembled beneath locations with 10–15 stalactites actively dripping water (see Supplementary Materials, Figure S2). PVC piping overlain with a plastic tarp makes up the frame, with a funnel attached to direct water to tipping bucket rain gauges (Onset HOBO RGB-M002) attached to micro-station data loggers (Onset HOBO H21-002) (Supplementary Materials, Figure S3). Due to the high amount of moisture in the cave passage, corrosion of rain gauge sensors led to instrument failure prior to 2012. This issue was addressed by retrofitting of the rain gauges with reed switches and pulse adapters in July 2012 [41]. The monitoring equipment was removed from the cave in September 2019.

2.2. Cave Drip and Precipitation Data

The cave drip dataset from James Cave consists of drip readings collected every 10 min from sites MS1 and MS2 (starting September 2007) and site MS3 (starting February 2008) until September 2019. Continuous time series data were offloaded from data loggers on a monthly to bimonthly basis until 2012, and bimonthly to seasonally until the instrumentation was removed in 2019. The data were processed using Aquarius Workstation [44] until 2014 and R Statistical Software [45] after 2014. Drip rates (mm per 10 min) were converted to discharge (mL/min) using a conversion of 3.7 mL to 0.2 mm of precipitation. In total, the drip dataset contains 427,486 observations for MS1, 459,248 observations for MS2 and 434,839 observations for MS3 [46]. Precipitation data were collected from the Climate Data Online (CDO) database [47]. The New River Valley Airport (Dublin, VA, USA) is the closest weather station to the entrance to James Cave (7.7 miles). The station location and the length of the precipitation record (1968–2021) were the factors which influenced the choice of precipitation data source.

2.3. Recession Methods

Recession analysis is traditionally used to distinguish between baseflow and surface runoff following a storm event, which is reflective of the hydrogeological properties in a watershed. In this study, we applied recession analysis on the drip hydrograph dataset to evaluate hydrologic characteristics of the epikarst and calculate storage volume in the epikarst. The recession equation, first derived by Boussinesq, Maillet and others [20] to examine stream discharge, is expressed as:

$${\mathrm{Q}}_{\mathrm{t}}={\mathrm{Q}}_0{\mathrm{e} }^{-\mathsf{\alpha }\mathrm{t}}$$

(1)

where Q₀ = the initial discharge (L³T⁻¹), Q_t = discharge at time t (L³T⁻¹), α = the recession constant (T⁻¹), and t = time (T). Plots of log Q over time can be used to distinguish different components (i.e., surface runoff and baseflow) of stream discharge hydrographs [48]. These shifts in slope, which reflect changes in α, have been used to interpret changes in storage or permeability. Initially, steep slopes reflect rapid flow; as the recession progresses, the slopes often become shallower, reflecting slower flow [28].

In addition to reflecting a transition between different flows or flowpaths, the recession constant α can also have a physical meaning, represented by:

$$\mathsf{\alpha }=\frac{\mathrm{K}}{{\mathrm{SL}}_{\mathrm{c}}}$$

(2)

where K is the hydraulic conductivity of an aquifer with length scale (L_c) and storage coefficient (S) [49]. Assuming that the storage and length scale of an aquifer remains constant, an increase in α should be associated with an increase in K. Conversely, as α decreases, K should decrease. Thus, changes in α should reflect a physical switching from flow paths of varying K [49].

There are other ways to express the recession equation, including [50]:

$$\mathrm{Q}={\mathrm{Q}}_0{\mathrm{k} }^{\mathrm{t}}$$

(3)

where k is the recession coefficient, defined as:

$$\mathrm{k}={\mathrm{e}}^{-\mathsf{\alpha }}={\left(\frac{\mathrm{Q}}{{\mathrm{Q}}_0}\right)}^{\frac{1}{\mathrm{t}}}$$

(4)

We used two approaches for conducting recession analysis of the cave drip datasets: a correlation method and a matching strip method. These methods, based on the recession equations outlined above, have been used previously to analyze recession data for spring discharge [24,26] and cave drips [5]. The correlation method, based on Equation (4), uses a plot of the current discharge against the discharge at a fixed lag time. This method plots an enveloping line of the correlated discharge observations where the slope defines the recession coefficient (k) for different lag times, t [28]. For each drip site, the correlation curve is the output with the slope of the envelope line defined as the recession coefficient k for a lag time of 1 day. Choosing the optimal lag time depends on the length of the recession periods. Previous studies [28] suggest selecting as long a period as possible that allows for a recession period to be analyzed. Because this method is focused on baseflow, we used this approach on subsets of the hydrographs with drip discharge below 50 mL/min reflecting drip baseflow. The optimal lag time was modeled by least squares regression to best fit the enveloping line for the curve.

The matching strip method involves plotting a full dataset of recession events in terms of time relative to the start of each event, and then fitting the data with the recession equation [51]. Recession events are ranked and sorted by the magnitude of decrease, after which point the dates are converted into relative time measurements in minutes from the beginning of recession. The segments (or strips) are adjusted horizontally in relative time to create a master recession curve (MRC), using a VBA program in Microsoft Excel [51]. This spreadsheet outputs an MRC for the drip record at all drip sites, which can be used to solve for the recession constant, α [34].

In addition to analyzing the larger datasets, we also analyzed multiple segments on individual drip hydrographs (hydrograph separation method), based on Equation (1) and following the approach from previous studies [24,28]. For this analysis, a select number of complete hydrographs were chosen for each drip site that had a hydrograph maximum drip rate >50 mL/min and a minimum drip rate of <10 mL/min. For the selected hydrographs, semi-log plots of Q vs. time were generated, and the slope of the line, reflecting α, was determined. If there was more than one slope identified in the semi-log hydrograph, the hydrograph was separated into two segments, generating an α for each slope. Each hydrograph was separated into two segments where there was a visible change in slope [24]. Those individual sections were then plotted as the slope of the natural log of discharge versus relative time. The α₁ and α₂ values were also used to calculate storage volume using:

$$\mathrm{V}\left(\mathrm{t}\right)=\underset{0}{\overset{\mathrm{t}}{{\displaystyle \int }}}\mathrm{Q}\left(\mathrm{t}\right)\mathrm{dt}=\frac{{\mathrm{Q}}_{0,1}}{{\mathsf{\alpha }}_1}\left(1-{\mathrm{e}}^{{\mathsf{\alpha }}_1\mathrm{t}}\right)+\frac{{\mathrm{Q}}_{0,2}}{{\mathsf{\alpha }}_2}\left(1-{\mathrm{e}}^{{\mathsf{\alpha }}_2\mathrm{t}}\right)$$

(5)

where V = storage volume of the epikarst draining to the cave drip (L³), Q_0,1 is the initial flow of the first segment of the recession curve (L³T⁻¹), α₁ is the recession constant of the first segment (T⁻¹), Q_0,2 is the initial flow of the second segment of the recession curve (L³T⁻¹), α₂ is the recession constant of the second segment (T⁻¹), and t is time (T) [24]. For this study, we calculated V₁ from the first hydrograph segment to reflect rapid flow from the epikarst and V₂ from the second hydrograph segment to reflect diffuse flow from the epikarst, as suggested in previous studies [24,28]. Using these values of V₁ and V₂, we then calculated the proportion of rapid flow as V₁/(V₁ + V₂) and the proportion of diffuse flow for the total volume for a drip event as V₂/(V₁ + V₂).

3. Results

3.1. Cave Drip Rate and Precipitation Records

The cave drip discharge data (2008 to 2018) collected in James Cave (Figure 2), along with the precipitation record for the same period (Figure 2; also see summary in

Table 1. Annual precipitation [47] for New River Valley Airport, Dublin, VA, USA.

Year	Precip (mm)
2008	537.7
2009	674.1
2010	619.8
2011	587.8
2012	389.6
2013	399.03
2014	1257.3
2015	1593.3
2016	1236.2
2017	1217.4
2018	1441.7

), show seasonal trends that are generally consistent from year to year. The seasonal trends include a period of increased drips (recharge), decreasing drips (recession) and no drips (dry period). The recharge period is defined as a period with high cave drip rates (>10 mL/min) and a marked response of cave drips to precipitation events [41]. The recharge period at the three drip sites typically starts in the early winter (January/February) but can start as late as March/April. In mid to late spring, the drip discharge goes into a period of recession, typically from April to August, where drips decrease exponentially. The dry period is defined as the period where drip discharge goes below 10 mL/min. This typically occurs in fall and winter, at which point the recharge season starts again.

3.2. Recession Analysis

3.2.1. Correlation Method

Figure 3 shows the results from the correlation of drip measurements to the previous drip measurement (lag drip rate). The slope of the envelope line fitted to the correlation curve is the recession coefficient, k, which was converted into α using Equation (4). Fitting of the envelope line was conducted to minimize residuals between the drip curves and the line. Linear regression analysis results are shown in

Table 2. Correlation method results for baseflow of MS1, MS2, and MS3, including recession coefficient (k) and R² for the best fit on the envelope. Conversion from k to α (1/day) was done using Equation (4).

Drip Site	k	Calculated α (1/day)	Calculated α (1/min)	R2
MS1	0.925	7.80 × 10−2	5.41 × 10−5	0.92
MS2	0.90	1.05 × 10−1	7.32 × 10−5	0.96
MS3	0.92	8.34 × 10−2	5.79 × 10−5	0.98

.

3.2.2. Matching Strip Method

Figure 4 shows the results for the MRC generated using the matching strip method (Excel spreadsheet provided by [51]). The fits of the exponential equation to the dataset have overall high R² values; however, there are deviations of the data from the exponential equation, especially at the higher flows for MS2. Summary results of the analysis are shown in

Table 3. Master recession curve results using the matching strip method for optimal fit over full drip record (2008–2018), Q₀ = discharge at the beginning of recession.

Drip Site	Q0 (mL/min)	α (1/Day)	α (1/min)	R2
MS1	147.26	1.73 × 100	1.20 × 10−3	0.99
MS2	223.67	7.05 × 10−1	4.88 × 10−4	0.99
MS3	702.8	3.41 × 100	2.37 × 10−3	0.98

.

3.2.3. Hydrograph Separation Method

Figure 5 shows results of the hydrograph separation method to determine α₁, and α₂ for individual hydrographs for each drip site. Note that the fits for α₁ are generally better than for α₂; at lower flows, there is lower accuracy of drip discharge measurement. Using the values of α₁ and α₂, the storage volume was calculated using Equation (5) for each segment (

Table 4. Storage volumes (L) and α values calculated for drip hydrographs for MS1, MS2 and MS3 using Equation (5). V₁ and a₁ reflect diffuse flow from epikarst; V₂ and α₁ reflect diffuse flow from epikarst. Also shown is the % rapid (concentrated) flow and % diffuse (base) flow for each hydrograph, using Equation (5) and the assumptions outlined in the text. Grey shaded regions reflect periods of recharge.

Drip Site	Date	V1 (L)	V2 (L)	V1 + V2 (L)	α1 (min−1)	α2 (min−1)	% Rapid Flow	% Diffuse Flow
MS1	1–28 April 2008 (recharge)	364.89	88.6	453.49	2.00 × 10−4	6.21 × 10−5	80.5	19.5
	8 July–30 August 2013 (recession)	207.5	1013.88	1221.38	7.71 × 10−5	3.50 × 10−5	17.0	83.0
	15 May–15 June 2018 (recession)	74.15	106.9	181.06	2.30 × 10−4	4.34 × 10−5	41.0	59.0
MS2	1–28 April 2008 (recharge)	163.96	36.46	200.42	6.30 × 10−4	5.69 × 10−5	81.8	18.2
	1–20 April 2009 (recharge)	837.63	198.31	1035.94	1.00 × 10−4	5.27 × 10−5	80.9	19.1
	1 July–10 August 2013 (recession)	728.44	419.3	1147.74	2.00 × 10−4	6.00 × 10−5	63.5	36.5
	25 April–1 June 2014 (recession)	868.89	208.32	1077.21	1.00 × 10−4	2.21 × 10−5	80.7	19.3
	20 April–15 June 2015 (recession)	485.87	182.35	668.21	1.20 × 10−4	2.45 × 10−5	72.7	27.3
	25 May–1 July 2018 (recession)	82.69	10.14	92.83	2.00 × 10−4	6.19 × 10−6	89.1	10.9
MS3	1–25 May 2009 (recharge)	369.18	19.52	388.69	1.70 × 10−3	3.79 × 10−5	95.0	5.0
	25 April–1 June 2014 (recession)	168.1	75.02	243.12	1.20 × 10−3	1.30 × 10−4	84.0	16.0
	20 May–1 June 2017 (recession)	105.51	6.73	112.24	5.00 × 10−4	2.75 × 10−5	69.1	30.9

). α₁ estimates were on average 3.64 times larger than α₂ for drip site MS1, for MS2 α₁ is on average 6.13 times larger than α₂, and for MS3, α₁ is on average 18.5 times larger than α₂. The variation of α₁ and α₂ is greatest for drip site MS3 and lowest for drip site MS1.

The hydraulic conductivity and storage properties of the epikarst can be characterized by the differences observed between α₁ and α₂. From Equation (2), α represents the ratio of hydraulic conductivity (K) and storage (S). In this study, two α parameters for rapid (α₁) and diffuse flow (α₂) were differentiated using the hydrograph separation method. As α₁ is greater than α₂, the ratio of hydraulic conductivity to storage is higher for α₁ than for α_2. Since storage is a physical property of the epikarst, during individual events, the average storage should remain constant during the event. For this reason, fluctuations in α likely reflect fluctuations in hydraulic conductivity of different portions of the epikarst (i.e., matrix vs. fractures) as the epikarst fills and drains.

4. Discussion

4.1. Diffuse vs. Rapid Recharge in James Cave

The three methods of recession analysis applied in this study can be used to address different types of flow (diffuse vs. rapid) contributing to the overall recharge from cave drips. The correlation method of hydrograph analysis focuses on baseflow [28]. We associate this baseflow with diffuse recharge, where flow is occurring mostly through pore spaces. The hydrograph separation method allows for separating components of the hydrograph into rapid flow (represented by α₁) and diffuse (represented by α₂). In contrast, the matching strip method used to develop the master recession curve results in average characteristics of the hydrograph, encompassing both rapid and diffuse components.

Combined, the recession analysis results (

Table 5. Values of α (min⁻¹) for MS1, MS2, and MS3, using the matching strip method (average flow), the correlation method (base/diffuse flow) and the individual hydrograph separation method (rapid flow and base/diffuse flow).

Drip Site	Matching Strip (Average)	Individual (α1) (Rapid)	Individual (α2) (Base/Diffuse)	Correlation (Base/Diffuse)
MS1	1.20 × 10−3	1.70 × 10−4	4.67 × 10−5	5.41 × 10−5
MS2	4.88 × 10−4	2.25 × 10−4	3.67 × 10−5	7.32 × 10−5
MS3	2.37 × 10−3	1.27 × 10−3	4.67 × 10−5	5.79 × 10−5

) show that the three drip sites have overall similar hydrologic characteristics, as reflected by similarities in the average α (10⁻³ to 10⁻⁴ 1/min) as derived from the matching strip method. Rapid flow α values (1/min) are in the 10⁻³ to 10⁻⁴ range, while diffuse are 1 to 2 orders of magnitude lower (10⁻⁵ range). At the broad scale, these values make sense in relation to one another.

The relative proportion of rapid flow to diffuse flow (see Table 4) suggests that during the recharge period, the dominant flow is rapid (>80%). For example, during the recharge period of 2009 at MS3 (1–25 May 2009) the proportion of rapid flow is 95% of the total with base flow making up 5%. In addition, an event during the recharge period of 2009 at MS2 (1–20 April 2009) showed 81% of flow was rapid flow with 19% as base flow.

In contrast, during drip recession, the dominant flow shifts to diffuse flow. For example, during an event at MS1 early in the recession period of 2018 (15 May–15 June 2018), the relative proportion of rapid to base flow is 41% to 59%. Another event in the recession period of MS2 in 2013 (1 July–10 August 2013) shows a relative proportion of rapid flow to base flow as 63.5% and 36.5%, respectively. This diffuse flow is likely through pore spaces in the soil and weathered bedrock of the epikarst. Although there are slight differences in the correlation plots between the three drip sites, the resulting k (and thus, α) values for diffuse flow are similar.

Comparing the α values estimated using the three methods shows that the matching strip method, which produces an overall α for many hydrographs, yields similar values (4.88 × 10⁻⁴ to 2.37 × 10⁻³ min⁻¹) to the a₁ values calculated using individual hydrograph analysis (1.7 × 10⁻⁴ to 1.27 × 10⁻³ min⁻¹) (Table 5), suggesting that rapid flow dominates the hydrograph. These values of a₁ are up to two orders of magnitude greater than the a₂ values calculated for diffuse flow using the correlation method (5.41 × 10⁻⁵ to 7.32 × 10⁻⁵ min⁻¹) and the individual hydrograph analysis (3.67 × 10⁻⁵ to 4.67 × 10⁻⁵ min⁻¹) (Table 5).

Each of the recession analysis methods utilized in this study has advantages and disadvantages. The generation of a master recession curve using the matching strip method is beneficial for evaluating the entire dataset, and its results reveal overall characteristics of the hydrologic system feeding the cave drips. The disadvantages are that the exponential model does not necessarily fit the entire dataset (see MS2, Figure 4) and that it does not account for variability that exists within the individual hydrographs. The correlation method is advantageous because the entire dataset can be analyzed, but the method is used explicitly to analyze base (diffuse) flow. The advantage of using the separation of hydrographs (α₁ and α₂) is that the method allows for examination of different flowpaths for individual hydrographs. Disadvantages of the separation method include that (1) the parameter estimates are for single hydrographs and not the full dataset and that (2) the hydrographs often show more than one change in slope, necessitating additional separation.

4.2. Conceptual Model of Storage and Recharge in the Epikarst at James Cave

The results from recession analysis and the variability characteristics of the cave drip discharge at James Cave from 2008 to 2018 suggest that the epikarst above James Cave contains a network of recharge pathways and storage zones that occur in pore spaces and fractures, supporting previous work at the site [7]. During recharge events in the winter and early spring, precipitation infiltrates the soil and weathered bedrock and starts to fill fractures and subsequently pore spaces. When the pore spaces connect and fractures “wet up”, infiltration can flow through these pathways to create the drip hydrographs that we see in the recharge period.

Drip sites MS1 and MS2 are similar in terms of the α values from recession analysis and variability characteristics, but there are some subtle differences. Comparison of storage volumes during the same recession event (1 July–10 August 2013) shows that MS1 had a higher proportion of diffuse flow (83%) than MS2 (37%), suggesting differences in hydraulic properties or storage characteristics of the pathways feeding the drips at MS1 and MS2 (Table 4).

Drip site MS3 is different from MS1 and MS2, with higher α₁ values (see Table 4 and Table 5). During drip events, the percentage of flow at MS3 is mostly rapid (up to 95%, see Table 4) with a lower proportion of baseflow. This can also be seen visually in the shape of the hydrographs with higher flashiness of discharge relative to MS1 and MS2.

Overall, these results suggest that cave drips at MS1 and MS2 have recharge contributions from both fracture and pore flow that originate from similar sources. In contrast, the recession analyses reveal that MS3 has a differing balance between rapid and diffuse flow than MS1 and MS2. In addition, the higher values of α suggest that concentrated flow through fractures is dominant in the epikarst that drains to MS3, which may reflect a more weathered or “karstified” epikarst. The differences in characteristics that we observe between MS1, 2 and MS3 also suggest heterogeneity in the permeability and storage characteristics of the soil and weathered bedrock in the epikarst of James Cave.

4.3. Study Limitations

The dataset used for this analysis has several limitations. First, the datasets are incomplete, as the instruments (rain gauge and datalogger) failed multiple times during the 10 years of data collection. The environmental conditions in the cave passage are marked by high humidity; thus, the rain gauges and data loggers used to collect the drip measurements were under constant exposure to conditions leading to corrosion. Due to the frequency of site visits in the early part of the study [39,41], the equipment was serviced or replaced to resume data collection within days of the instrument failures. However, after 2013, there were fewer site visits and thus instrument failures were often not recognized for several months. In addition, there are periods of missing data during high drip events, when the drip tarps would fill with silt, clogging the rain gauges. This prevented comparison of drip rates in periods of low precipitation (2008–2013; 400–700 mm annual precipitation) to periods of higher precipitation (2014–2018; >1000 mm annual precipitation) (see Table 1). Lastly, the study is limited to only three drip monitoring stations. The monitoring stations had to be carefully selected for this study to ensure ease of access for data downloads. Thus, inferences made in this study are limited to the data recorded individually at each of the three monitoring stations. Future research could address some of these data limitations by using more robust monitoring methods to measure drip rates or utilize interpolation methods to fill in the gaps in the time series.

There are also limitations in the methods used in this study. First, the recession analysis methods assume homogeneous and isotropic media, and Darcian (pore) flow. However, results from this study suggest that permeability varies throughout the epikarst. Thus, although we needed to make these assumptions to do the analysis, we know that the assumptions are likely violated. There are other equations that are used for fracture flow, but there are not equivalent recession analysis methods for these equations. Second, hydrograph separation methods to delineate distinct segments can also be challenging as it can be difficult to decide how to make the delineation. In this study, the delineation was done visually, which is subjective and has associated errors. Third, correlation of drip rate to lag rate is limited in that it assumes linearity between subsequent calculations. This is a necessary assumption to fit an envelope line to the curve which determines the recession coefficient [28]. Last, using α to evaluate hydrogeologic characteristics is imperfect, as it assumes that there is a physical meaning to this parameter. As shown in Equation (2), α is a lumped parameter that reflects the ratio of hydraulic conductivity (K) to storage coefficient (S). Thus, using α to interpret K or S individually can be challenging. Future research could address some of the method limitations by utilizing equations specifically for fracture flow.

5. Conclusions

This study utilized a 10-year cave drip record to evaluate recharge and storage in the epikarst of a local cave (James Cave, Pulaski, VA, USA). Three recession analysis methods were utilized: a correlation method, a matching strip method, and a hydrograph recession method. Results of the recession analysis show a cave drip response that reflects both rapid recharge through fractures and diffuse flow due to draining of pores. Results show a similarity in hydrologic characteristics, reflected by the recession constant α and drip variability characteristics, at two of the drip sites located within 100 m of each other (MS1 and MS2), which are different from a third site (MS3). Overall, this study suggests that temporal and spatial characterization of cave drip rates is important for quantifying recharge to karst aquifers.