The Edwards Aquifer (EA) in Central Texas provides water supply for over two million people and contains springs that are hydrologically and ecologically important to the region. The residence time of groundwater in the EA ranges from a few days to many thousands of years, since water in the aquifer is contained and transported within both matrix porosity and large conduits. In this study, stable isotopes of water from five springs are investigated for tracing the origin of water and hydrological processes in the EA system during 2017–2019. There is a quick response of the isotopic signals measured at these springs to changes in the isotopic compositions of precipitation. By utilizing an isotope mixing model, we have identified sources of water for these springs with a bi-modal distribution of groundwater supply in the EA: water supplied from deep groundwater with a longer residence time (an average of 67%) and supplemental epikarst interflow with a shorter residence time (an average of 33%). The evolution of hydrochemical water types from HCO3–Ca to HCO3·Cl–Ca·Mg along the EA flowpaths indicates that inputs from epikarst interflow are greater in springs within the artesian zone than the springs within the contributing zone.
Stable isotope tracers were utilized to understand the origins of spring water and hydrological processes in the Edwards Aquifer.
There is a quick response for isotopic signals in spring water to changes in the isotopic compositions of precipitation.
There is elevated input from epikarst interflow for spring water balance in the artesian zone than the contributing zone.
The Edwards Aquifer (EA), located along the southeastern edge of the Edwards Plateau in Texas, is one of the most prolific artesian aquifers in the world (Sharp & Banner 1997). The EA is the major source of potable water for two million people and the primary water supply for agricultural and industrial use in Central Texas (Schindel et al. 2004). The EA is primarily recharged by surface runoff that feeds the numerous streams and creeks originating in the Texas Hill Country. The primary natural discharge points from the EA are springs in the artesian zone, especially those along the Balcones Fault Zone (Hunt et al. 2019). These springs have supported natural ecosystems, human populations, and urban development in Central Texas for thousands of years. However, water shortages can occur in the EA due to excessive pumping of groundwater and periodic droughts (Tian et al. 2020).
The EA confined/artesian zone contains highly variable (in time and location) mixtures of waters with residence times ranging from a few days to many thousands of years (Eberts et al. 2011). Radiocarbon dating is not a reliable groundwater dating method in karst aquifers (Maloszewski & Zuber 1991) because the contribution of ‘dead’ carbon from the dissolution of calcite/dolomite would lead to a significant over-estimation of the groundwater's age. Hunt et al. (2016) determined apparent ages for well water from the EA by applying a 3H/3He dating method, which indicated that the groundwater is the result of mixing of both modern (less than 60 years) and premodern (higher than 60 years). Kuniansky et al. (2001) utilized the finite-element method to model travel times for water within the aquifer and estimated a range of 350–4,300 years, assuming a flowpath from the West Nueces River Basin to Comal Springs. Dye tests in the EA usually reveal high velocities and short residence times. For example, Hunt et al. (2005) utilized fluorescein dye to determine that the water arrived at Barton Springs between 7 and 8 days from the recharge cave.
The rocks hosting the EA contain large voids such as conduits, caves, fractures, and fissures, as well as matrix porosity (Hovorka et al. 1995; Chen et al. 2012), which facilitate recharge of the aquifer due to precipitation. Hydrographs from wells and/or springs can be used to analyze the structure and response of the hydrogeological system (Halihan & Wicks 1998). However, the range of mechanisms affecting hydrographic responses has not been thoroughly investigated in the EA region (Zhang et al. 2020). Hydrochemical and isotopic analysis of spring water in a watershed are particularly suitable for tracing the origin of water and providing information on hydrological processes (e.g. Araguás-Araguás et al. 2000; Marfia et al. 2004; Martinez et al. 2015; Shakya et al. 2019). This research has investigated isotopic signatures of five springs from the EA region through monthly sampling during 2017–2019. Changes in the isotopic compositions of the springs were compared with the isotopic signatures in local precipitation. The objectives of this study were: (1) to document the spatial and temporal variations in the isotopic composition of spring water and precipitation, (2) to investigate the isotopic response of spring water to the isotopic variations in precipitation, and (3) to determine the water mass balance for different sources potentially feeding springflow.
The EA is an extensively karstified groundwater flow system developed in limestone rocks deposited during the Cretaceous (Schindel 2019). As shown in Figure 1, the EA system has been geographically divided into three major hydrogeologic zones: the contributing zone (or drainage area), the recharge zone, and the artesian zone (Stone & Schindel 2002). The contributing zone collects rainwater and springflow, which forms the headwater of major streams in the region; then the surface runoff is conveyed into the recharge zone without significant loss. The recharge zone approximately overlies the Balcones fault zone, allowing for vertical downward movement of surface water into the EA. The artesian zone, where the EA is confined, contains most of the large production wells that pump water out of the EA. The general flowpaths of water within the EA are from the north to the south in the contributing zone, then from the southwest to northeast in the artesian zone (Figure 1(a)).
The hydrogeologic cross section within Bexar and Comal Counties is shown in Figure 1(b). The Blue Hole (BH) is the largest spring of ‘The San Antonio Springs’, which is located at a nature preserve called the Headwater Sanctuary on the campus of the University of Incarnate Word. The San Pedro Springs (SPS) is located a few miles northwest of downtown San Antonio, directly adjacent to the campus of San Antonio College. The Comal Spring (CS) is the largest spring in Texas, and it is located at the base of a steep limestone bluff in Landa Park of New Braunfels. The altitudes for BH, SPS, and CS are 206, 202, and 192 m, respectively. The J-17 Index Well is located at Fort Sam Houston in San Antonio. The J-17 Well is on a major flowpath in the EA and responds quickly to pumpage and recharge. Thus, the water level variation of J-17 Well is a general indicator of the amount of water within the EA. Springflow at CS becomes intermittent when the water level measured at J-17 Well drops below 192 m; flow at BH ceases when the J-17 Well level is below 207 m; and flow at SPS ceases when the J-17 Well level is below 201 m. No flow was recorded at BH and SPS several times during our sampling interval. Water was constantly issuing from CS during our sampling interval; however, historically, flow ceased during the 1950s drought when the J-17 Well dropped to 187 m (Cox et al. 2009).
Climatic and weather conditions
The research area, located in Central Texas, has a transitional humid subtropical climate, with hot summers and generally mild winters with infrequent severe freezes. The climatic parameters (temperature and precipitation amount) at the City of Boerne during 2017–2019, together with those from 1981 to 2010 for comparison, are summarized in Figure 2. The weather station site is within 20 km from our precipitation sampling site, and the precipitation falling on the contributing zone should be the source of water for the springs in the artesian zone.
There were high summer temperatures and increasing seasonal variabilities of precipitation during 2017–2019, compared with long-term averages in the period 1981–2010 (Figure 2). Monthly precipitation amounts varied between 53 mm (January) and 118 mm (May) derived from the long-term (1981–2010) monthly means, and the monthly temperature varied between 9.0 °C (January) and 27.6 °C (August). Notably, the precipitation amounts in the summer seasons (June–August) were significantly lower than the corresponding averages. However, September 2018 was the wettest September on record in Boerne: 360 mm in September 2018 compared with an average of 88 mm during 1981–2010. The average annual precipitation amount was 964 mm during 1981–2010. Annual precipitation amounts for 2017 (733 mm) and 2019 (672 mm) were lower than the average, whereas the annual precipitation amount for 2018 (1,092 mm) was higher than the average. The annual mean temperatures during 2017–2019 were 20.3, 19.4, and 19.3 °C respectively, which were higher than an average of 18.8 °C from 1981 to 2010.
The water levels of the J-17 Well, discharge of CS, along with meteorological data during 2017–2019 are shown in Figure 3. The springflow of CS originates from groundwater flow in the southwest (Figure 1), and the flowpath moves past the J-17 Well on its way toward CS. The underground flow might take a few years based on MODFLOW calculations for the deep groundwater with a long residence time. However, the hydrograph illustrates both J-17 Well levels and CS springflows quickly respond to net precipitation (precipitation–evaporation), and there is almost no lag between CS springflow and J-17 Well levels (Figure 3). The fast response to net precipitation for both sites should be due to additional inputs from surface runoff or epikarst interflow. Our spring water are collected at the spring orifices to avoid inputs from surface runoff or stormwater. There are typical karst features of fissure, fracture, sinkholes, and caves near the spring sites, which facilitate surface water infiltration. Therefore, we conclude that the additional supplement for the spring water balance is epikarst interflow or shallow groundwater with short residence time.
Water sampling and field measurements
This research has investigated isotopic signatures of five springs from the EA, two springs (BS1 and BS2) from the unconfined section and three springs (BH, SPS, and CS) from the confined section, by monthly sampling during 2017–2019, along with the isotopic signatures in precipitation. Locations, altitudes, and site descriptions of the five investigated springs are listed in Supplementary Material, Table S1. Springs termed BS1 and BS2 emerge from the contributing zone and represent the headwater of the Upper Cibolo Creek (UCC) watershed. The UCC generally flows in a southerly or easterly direction across the UCC watershed. Two major springs, BH and SPS, discharge from the artesian zone proximal to the downtown region of San Antonio. CS lies within the transition of the recharge zone/artesian zone. Precipitation samples were collected after each rain event in an open area near the main campus of UTSA (29°35′N, 98°38′W). Locations and photos for these sampling sites are shown in Figure 4 and Supplementary Material, Figure S1, respectively.
A total of 148 spring samples were collected during 2017–2019: 36 samples for both BS1 and BS2, 32 samples for CS, 29 samples for SPS, and 15 samples for BH. All investigated springs were sampled monthly in general; however, spring discharge of BH and SPS had ceased many times during our sampling interval.
Samples collected for stable isotope, cation, and anion analysis were filtered using 0.2-μm filters at the time of collection and stored in sterilized HDPE bottles. The bottles and caps were rinsed three times with filtered water, sealed with parafilm, and refrigerated until the time of analysis.
Water temperature, electrical conductivity (EC), total dissolved solids (TDS), and pH were measured at the time of collection. Water temperature, EC, and TDS were measured with a YSI® 556 MPS handheld multi-parameter sensor, while pH was measured with a Hach® sensION1 handheld pH meter. EC and pH sensors were calibrated before each sampling field.
Daily rainwater was collected in a 10-cm diameter high capacity Stratus® rain gauge. A ping-pong ball was placed at the mouth of the rain gauge to prevent evaporation (Michelsen et al. 2018). Sample collection was usually performed during the morning at about 8 a.m. when rain events happened at nights or early mornings, while some collections were performed at other times when the rainfall occurred during the daytime. The rain gauge was emptied and cleaned after each rain event, and sub-samples of the homogenized precipitation were sealed in pre-cleaned bottles, with no headspace, for future δ18O and δD analysis (Gröning et al. 2012).
The concentrations of cations (Ca2+, Mg2+, Na+, K+, and Sr2+) and anions (, Cl−, F−, Br−, and ) were measured using Dionex ion chromatography (IC) at Texas State University. Water samples with conductivities higher than 0.90 mS/cm were diluted with Milli-Q ultrapure water in a specific ratio before analysis. Milli-Q blanks and standards with dilutions bracketing the range of ion concentrations in the samples were run at the time of analysis.
The alkalinity of filtered water samples was measured by titration with 0.01 mol L−1 hydrochloric acid. The pH of our collected groundwater samples generally ranges between 6 and 9, which implies the concentration is negligible, and the alkalinity should be approximately the same as concentration. Thus, we only tested the bicarbonate alkalinity, i.e. concentration, using the bromocresol green-methyl red mixed indicator.
Isotope mixing model
Based on the analysis of groundwater hydrographs (Figure 3), we speculated that there are both shallow and deep groundwater contributions in the spring discharge mass balance (Hare et al. 2021). In the EA context, the shallow groundwater has shorter residence time which could range from days to months for ‘quick’ response; the deep groundwater has longer residence time which could range from years to centuries for ‘slow’ response (Hunt et al. 2016).
RESULTS AND DISCUSSION
Hydrochemical data visualization
The field measurements, major ion concentrations, and isotopic analysis for all samples of five investigated springs during 2017–2019 are listed in Supplementary Material, Table S2. The summary statistics of hydrochemical and isotopic parameters of the five springs are listed in Table 1, for each spring is listed separately in Supplementary Material, Table S3. We utilized major cations and major anions to calculate the charge balance error (CBE), and the CBE is mostly clustered within ±5%, with many outside of ±5% and a few outside of ±10%. Thus, the hydrochemical dataset is acceptable for general hydrochemical characterization.
The temperature and pH are consistent among the five springs, with the lowest coefficient of variation (CV), indicating that the major source of the spring water is from the deep underground aquifer. Among water isotopes, D-excess (29%) is more variable than δ18O (15%) and δD (11%). The range of concentrations for Ca2+ and are relatively small, with relatively low CVs for Ca2+ (28%) and (33%). The concentration of these ions results from the dissolution of limestone within the EA (e.g. Engesgaard & Christensen 1988; Tian et al. 2020), and the spring water are likely supersaturated generally with carbonate in the study area. The concentrations for K+ and display substantial variabilities among the spring sites, with relatively high CVs for K+ (142%) and (49%). None of these ions is from the dissolution of limestone in the research area but is more likely associated with human activities (e.g. Jiang et al. 2008; Wang et al. 2018). Their highly variable concentrations indicate the sources of water not only come from deep groundwater, but with additional supplementary inputs from precipitation, surface runoff, or epikarst interflow.
The correlation coefficients of chemical parameters for the spring samples are shown in Figure 5. There are significant correlations among ion concentrations of , Cl−, , Mg2+, K+, and Na+. Ion concentrations of Mg2+ and likely originated from the dissolution of gypsum, while the other ions are related to human activities. Previous studies have indicated elevated chlorine and nitrate levels for surface water and groundwater in the research area (e.g. Morrall et al. 2004; Sullivan & Gao 2016). A SWAT model coupled with groundwater flow and contaminant transport models indicated that anthropogenic inputs as well as biogenic wastes were important sources of nitrate to the EA system (Sullivan et al. 2019). D-excess is negatively correlated with δ18O and δD. The isotopic variations are well correlated with ion concentrations of and Cl−. The possible reason is that both the isotopic signatures and contaminants of these springs are precipitation input signals. The higher and Cl− concentrations are not directly from precipitation inputs, but the load of anthropogenic contaminants during the precipitation infiltration.
A Piper diagram displaying relative concentrations of cations and anions for the five springs is shown in Figure 6, and Piper diagrams for each spring are shown separately in Supplementary Figure S2. The relative concentrations of anions and cations for these spring samples are shown in the Schoeller diagram (Supplementary Material, Figure S3). The dominant species in the spring water are Ca2+ and for all five springs, which indicates that the major contribution of ions originates from water–rock interaction within the aquifer. Thus, the water source of these springs is mainly from ‘older’ deep groundwater. However, occasionally there are elevated concentrations of Cl−, which may indicate additional supplementary inputs from ‘young’ water, more directly tied to epikarst interflow or shallow groundwater (Moore et al. 2009; Chaudhuri & Ale 2014; Retike et al. 2016). There are higher Cl− concentrations in the three springs (BH, SPS, and CS) from the confined zone than the two springs (BS1 and BS2) from the unconfined zone. The hydrochemical classification also evolves from HCO3–Ca type to HCO3·Cl–Ca·Mg type from the contributing zone to the artesian zone, which indicates the springs in the artesian zone might have high contributions of precipitation input.
Isotopic spatiotemporal variations
The cross plots of δ18O and δD for the spring water and precipitation are provided in Figure 7; the cross plot only for the spring water is provided in Supplementary Material, Figure S4. The local meteoric water lines (LMWLs) are derived by the monthly and daily precipitation isotopic data using ordinary least squares regression (OLSR) (Crawford et al. 2014). The groundwater line (GWL) fits the data from the springs and plots below the LMWL, which are commonly interpreted as the characteristic of residual waters from partial evaporation (Bowen et al. 2018). The slope and intercept of LMWLs are close to each other; thus, monthly isotopic variability still preserves the isotopic signals of rain events. Most spring water samples fall close to the LMWL, indicating that the spring discharge undergoes minimal evaporation before precipitation recharge in general. A few spring water samples, especially from the three springs within the artesian zone, plot to the right of the LMWL with a decreased slope, which indicates a certain extent of evaporation for these samples. All the spring water are located on the lower part of the monthly precipitation, which suggests that the groundwater from the EA generally has relatively lighter isotopic signatures than the precipitation signals during 2017–2019.
The global meteoric water line (GMWL) is defined as δD = 8δ18O + 10 (Craig 1961). The slopes and intercepts of LMWLs are slightly lower than the GMWL, which are controlled by local climatic factors. The intersection points between LMWLs and GWL are more negative compared with the amount weighted mean δ18O (−2.7‰) in precipitation during 2017–2019. Long-term mean δ18O in precipitation (from 1950 to present) in the study area is obtained from the Online Isotopes in Precipitation Calculator (OIPC, http://www.waterisotopes.org). The modeled long-term mean δ18O in precipitation (−4.3‰) is close to the intersected δ18O values between LMWLs and GWL, which indicates the theoretical initial spring recharge originated from the deep groundwater with longer residence time.
Violin plots of the isotopic variations for the springs and precipitation are shown in Figure 8. Tukey multiple comparisons for the means of isotopic signatures in the spring water and precipitation are provided in Supplementary Material, Figure S5. There is no significant difference for δ18O and δD among the three springs in the artesian zone (around −4.0‰ for δ18O), whereas two springs in the contributing zone have lighter isotopic signatures (around −4.5‰ for δ18O). The δ18O and δD signals in precipitation have generally heavier isotopic signatures than the spring water. Thus, the springs from the contributing zone and artesian zone might have different contributions of deep and shallow groundwater. A few outliers above the δ18O range for the springs might be the response to precipitation input or epikarst interflow. There is no significant difference for D-excess between the springs and precipitation (around 10‰). D-excess has been widely used to identify source regions of precipitation (e.g. Cui et al. 2009; Lewis et al. 2013; Steen-Larsen et al. 2015). However, D-excess might not be a good tracer for the springs in the mass balance calculations due to its significant variations. Thus, δ18O of these springs is utilized to examine the possible precipitation input for the spring water balance later.
The temporal variations of δ18O for the spring water and precipitation, together with the observed monthly precipitation amount in San Antonio, are shown in Figure 9. The amount weighted mean δ18O of precipitation for all rain events during 2017–2019 is −2.7‰. The δ18O of spring water shows small variations and varies below the mean δ18O of precipitation during most of the time during 2017–2019. This demonstrates the bi-modal distribution of groundwater age for these three springs: the springs mainly originate from deep groundwater, with additional supplemental inputs from epikarst interflow or shallow groundwater.
The modeled long-term mean δ18O in precipitation (from 1950 to present) is −4.3‰, and a weighted average δ18O of −4.1‰ in precipitation (1999–2007) was obtained in a nearby city Austin, Texas (Pape et al. 2010). The springs investigated in this work display δ18O signatures close to long-term mean δ18O in precipitation, which indicates that these springs mainly originated from the deep groundwater with longer residence time. However, the amount weighted mean δ18O for all rain events during 2017–2019 is −2.7‰, which is significantly higher than the mean δ18O values observed in our spring sites and previously observed δ18O in precipitation from Texas. The enriched δ18O values in precipitation during 2017–2019 might be due to higher annual mean temperatures or increasing seasonal variability of precipitation.
The isotopic signatures for the five investigated springs demonstrate similar trends; however, the springs in the contributing zone (BS1 and BS2) generally have lighter isotopic signatures than the springs in the artesian zone (CS, BH, and SPS). The isotopic signatures of precipitation have higher variations than the spring water. There is no lag or 1-month lag of this response for δ18O in springs to δ18O in precipitation. Thus, there is additional ‘young’ precipitation input for springs beside the ‘old’ deep groundwater, and this ‘young’ shallow groundwater could reach these spring sites within 1 month. The response of δ18O in springs to precipitation input is also related to the precipitation amount. For example, June and August 2018 have the heaviest isotopic signals, but there is no response of isotopic signals in June 2018 due to its small precipitation amount. The precipitation input is likely from the epikarst interflow, which agrees with its quick response and short-lasting time.
Mass balance calculations
Since all springs have consistent isotopic signals (also the lowest δ18O values) during their base flow conditions (Figure 9), we directly use the lightest isotopic values for the investigated springs as their endmember of δ18ODG for deep groundwater: −5.1‰ (BS1 and BS2); −4.6‰ (BH and CS); −4.4‰ (SPS). Since most spring water undergoes minimal evaporation before precipitation recharge (Figure 7), and ‘young’ shallow groundwater could reach these spring sites within 1 month (Figure 3), δ18OSGN for shallow groundwater could be directly assigned as δ18O in precipitation for each monthly interval calculated by Equation (4). Our observed isotopic data of daily precipitation are consistent with those observed in a previous event-based isotopic study of precipitation during 2015–2017 in Austin (Sun et al. 2019). We assumed that δ18O in precipitation might be spatiotemporally homogeneous in Central Texas and assigned the same δ18OSG for the investigated springs. In that sense, the altitude effect on the isotope in precipitation is considered limited in the research area. The altitude of the rain collector is 303 m, about 170 m lower than the two springs (BS1 and BS2) from the unconfined zone and 100 m higher than the three springs (BH, SPS, and CS) from the confined zone. However, the altitudinal lapse rate is unknown in the study area, which is the primary source of uncertainty in our isotope mixing model.
The monthly FSG results for the five springs are shown in Figure 10(a). Note that the isotope mixing model would not return FDG and FSG contribution estimates when the isotopic values for spring water and the endmembers have no significant differences. The amount weighted mean δ18O in precipitation during 2017–2019 is −2.7‰, which could be utilized as the endmember of mean δ18OSG during 2017–2019, assuming that all precipitation events lead to recharge. The overall FDG and FSG results for the five spring during 2017–2019 are shown in Figure 10(b). The monthly FSG for all springs display broad ranges of change from 0 to 80%. The monthly FSG for the five springs demonstrate similar trends during 2017–2019, which verified the isotope mixing model's validity. The overall FSG during 2017–2019 varied from 20 to 40% (an average of 33%); however, the three springs in the artesian zone (an average of 39%) have higher overall FSG than the two springs in the contributing zone (an average of 25%). The relatively high contributions of precipitation input for spring in the artesian zone might be due to their relatively lower elevation and larger catchment area than the springs in the contributing zone. This elevated input from epikarst interflow for springs within the artesian zone agrees with the evolved hydrochemical types from HCO3–Ca type to HCO3·Cl–Ca·Mg type from the contributing zone to the artesian zone, which verified the validity of the isotope mixing model.
This research has investigated isotopic signatures of five springs from the EA region, two springs from the unconfined section and three from the confined section, by monthly sampling during 2017–2019. There is a quick response for isotopic signals in spring water to changes in the isotopic compositions of precipitation. By utilizing the isotope mixing model, we have identified sources of water for these springs with a bi-modal distribution of groundwater age: mainly originate from deep groundwater (an average of 67%), with additional supplement shallow groundwater (an average of 33%). With the evolved hydrochemical type from HCO3–Ca type to HCO3·Cl–Ca·Mg type along the EA flowpaths, there is elevated input from epikarst interflow for spring water balance in the artesian zone than the contributing zone.
There is a lack of quantitative analysis on ages of the spring water due to no proper dating methods. Our isotope mixing model could untangle groundwater mixing ratios (shallow vs. deep groundwater) and travel times (short vs. long residence time) within the EA. However, there are several caveats in our isotope mixing models: the altitude effect on the isotope in precipitation is considered limited in the research area; we assume that spring water generally undergoes minimal evaporation; the monthly sampling of spring water may not be representative of temporal isotopic variations. More systematical investigations on the evolution of isotopic signatures from precipitation to groundwater are needed to refine the isotope mixing model in future research.
We thank Stuart Brown, Gary Rogers, Erica Martinez and Yuanxin Qu for their assistance with property access and water sampling. Gabrielle Timmins and Ashley Cottrell are thanked for their contributions to analyzing dissolved ions in the Texas State University Laboratory. The authors thank Dr Przemyslaw Wachniew (Associate Editor) and two anonymous reviewers for their valuable comments and suggestions. The activities reported here were funded in part by the Edwards Aquifer Authority (EAA) under Interlocal Cooperation Contract No. 19-928-AMS, but this work does not necessarily reflect the view or regulatory position of the EAA. This research was also funded in part by the National Key Research and Development Program of China (2016YFA0600504 and 2020YFA0607700), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB26000000), and the National Natural Science Foundation of China (41888101 and 41690114).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.