Low Flow Trends in Texas Stream Segments Serving Unique Hydrologic Functions

THE TEXAS WATER JOURNAL is an online, peer-reviewed, and indexed journal devoted to the timely consideration of Texas water resources management, research, and policy issues. The journal provides in-depth analysis of Texas water resources management and policies from a multidisciplinary perspective that integrates science, engineering, law, planning, and other disciplines. It also provides updates on key state legislation and policy changes by Texas administrative agencies. For more information on on the Texas Water Journal as well as our policies and submission guidelines, please visit texaswaterjournal.org. As a 501(c)(3) nonprofit organization, the Texas Water Journal needs your support to provide Texas with an open-accessed, peer-reviewed publication that focuses on Texas water. Please consider donating. The Texas Water Journal is published in cooperation with the Texas Water Resources Institute, part of Texas A&M AgriLife Research, the Texas A&M AgriLife Extension Service


INTRODUCTION
Low flows in streams are often used in natural resource management and environmental regulation as an indicator of overall stream health. Several researchers (e.g., Hisdal et al. 2004;Jowett and Biggs 2006;Bradford and Heinonen 2008;Thomas et al. 2019) have highlighted increased stress on aquatic, riparian, and hyporheic ecosystems during low flow due to decreased water availability and habitat quality. During these intervals, low flows help maintain longitudinal connectivity in the stream (Curran et al. 2012). Changes in flow and groundwater levels due to precipitation and seasonal factors have ecological impacts on stream communities. For instance, fish in riffle or shallow-water habitats can experience habitat loss. Low flows are critical for successful reproduction as many fish species migrate upstream to spawning sites (Bradford and Heinonen 2008;Bogan et al. 2017). In water resource management, Smakhtin (2001) noted that the evaluation of low flows is necessary for water allocations for competing interests such as municipal supply, irrigation, and recreation. The impact of low flow characteristics on water availability and water security has been discussed by Stahl et al. (2008), Vorosmarty et al. (2010), and Brauer et al. (2015). From an environmental health perspective, low flows such as the 7Q 10 (the annual minimum 7-day mean flow with a 10-year return period) have been used for prescribing total maximum daily loads (TMDLs) in some parts of the United States (Steinschneider and Brown 2012). In Texas, the 7Q 2 (7-day, 2-year low flow) is used to establish water quality criteria for wastewater discharges as part of the Texas Pollutant Discharge Elimination System (TPDES) by the Texas Commission on Environmental Quality (TCEQ;TCEQ 2010).
In recognition of the critical roles that streams play in local ecosystems, the 80th Texas Legislature passed bills (e.g., Senate Bill 3 2007) to develop, manage, and preserve the water resources of the State and protect instream and freshwater inflows. The 16 regional water planning groups in Texas recommended that the State Legislature designate certain segments of streams as ecologically significant. Streams with this designation are acknowledged for their unique ecological value in serving various biological and riparian functions, for supporting endangered or threatened species and communities, and for serving important hydrologic functions, including flow stabilization and groundwater recharge (Texas Administrative Code (TAC) § 358.2 2012). Upon receiving this designation, the stream segment is protected from the construction of State-funded reservoirs. Most of the segments that have been recognized for their unique hydrologic function are located in the Edwards Plateau (EP) region of Texas and serve as above-ground recharge sources for the Edwards Aquifer, one of the most productive karst aquifers in the world (Thomas et al. 2019). The majority of the Edwards Aquifer recharge (85%) is contributed locally by the overlying watershed while the remainder is sourced from direct precipitation and below-ground flows from adjacent aquifers (Edwards Aquifer Authority [EAA] 2003). A brief review of existing literature on low flows in Texas streams is presented to substantiate the motivation behind the current study. Poshtiri and Pal (2016) used various indicators to study the magnitude, timing, and duration of low flows in the continental United States from various starting periods to 2012. This study included several Hydro-Climatic Data Network (HCDN) stations in Texas. In general, they found a drying trend from 1980 onward (relative to pre-1980) for the Texas Gulf region. A few sites in South Texas also showed statistically significant decreasing trends in annual 7-day minimum flows. Without reference to any specific site or river basin, they reported an increase in the frequency of dry days in the Texas Gulf region. Thomas et al. (2019) used a diverse suite of hydrological indicators to assess hydroclimatic trends in ecologically significant stream segments in the Nueces River Basin from 1970 to 2014. They reported decreasing trends in annual minimum and annual median flows in four of the six gages used in the study with no corresponding conclusive trends in precipitation. They concluded that even small changes in land use and land cover in this basin, coupled with the lack of statutory oversight on water withdrawals in these segments, likely contributed to the declining trends in annual low flows. Recently, Rogers et al. (2020) evaluated trends in streamflow at selected locations in Texas as part of a larger study encompassing the southern and southeastern United States. They found statistically significant decreasing trends in flow for 1970-2015 at many 'reference' sites (i.e., sites with minimal anthropogenic influence) in Texas and concluded that these declines may be partially climate-driven. They also highlighted the year 1970 as being the beginning of a period of significant decline in mean streamflow and noted that analyses beginning with this year may be useful for studying climate change impacts on streamflow.
Climatic variability, specifically increasing surface water temperature (which is a function of ambient air temperature) and potential evapotranspiration (PET), has been demonstrated to influence the health of individual taxa and ecosystem functioning in streams, particularly the cycling of carbon. In situations of drying or drought, streams provide critical drought refuges such as (1) remnant or perennial pools and seeps for surface water habitats; (2) sediment or stones for resting stages, and (3) the hyporheic zone for taxa capable of vertical migration Bogan et al. 2017). These refuges can provide a critical source to support biodiversity down-stream if located at headwaters (Bogan et al. 2015). However, not all taxa benefit from refuges and site-specific factors, such as substrate type or oxygen concentration. In addition, fluctuations in the pH conditions of pools can result in community structures (e.g., surface invertebrates such as aquatic insects) that are significantly different before and after recovery from drought (Acuña et al. 2005;Bogan et al. 2015). Even though some surface invertebrates do not find refuge in the hyporheic zone, both surface insect fauna and hyporheic non-insect fauna demonstrate overlap between intermittent and perennial streams (Del Rosario and Resh 2000). This suggests the ability of these fauna to re-colonize streams that are hydrologically connected via swimming, crawling, or flying (Bogan et al. 2017).
It is evident from the preceding literature that much interest has been shown in investigating the impacts of climate (both past and future) on streamflow. When watersheds undergo changes in land use or experience anthropogenic modifications, we often suspect that flow regimes may be altered and studies investigating the impact of these changes are often carried out. However, where ecologically significant stream segments are concerned, hydro-meteorological changes in their watersheds often go unnoticed. As highlighted earlier, these segments are of critical importance to the ecosystem, and yet very little literature exists on long-term changes in their hydrology. Therefore, there is an urgent need for studies that examine climate impacts on streamflow in these segments, particularly low flows. Smakhtin (2001) emphasized the need for such studies to receive more focus and recommended the use of a variety of low flow indices to understand this "dynamic concept." Therefore, overarching goal of this study is to examine trends in low flows at seven United States Geological Survey (USGS) gaging stations in Texas using a suite of complementary hydroclimatic indicators. These stations were selected due to their location on stream segments that serve a unique hydrologic function. The selected gages are maintained by the USGS as part of the HCDN (Lins 2012) and are minimally impacted by anthropogenic factors. Specifically, we use metrics that reflect the magnitude, duration, and frequency of various types of low flows; examine the concurrent trends in associated variables such as temperature, precipitation, and evapotranspiration; and develop drought indices to evaluate low flow trends in those segments holistically and identify potential drought drivers. To the best of the authors' knowledge, there have been no prior efforts to characterize low flow trends in segments serving such valuable hydrologic functions in Texas. Therefore, the discussion of our results in the context of potential meteorological drivers and ecohydrological implications is a novel feature of the study.

Study area characteristics and data overview
The State of Texas comprises ten distinct climatic divisions (CDs; Figure 1; National Climate Data Center [NCDC] 2015). Regions that fall within the same CD are similar in seasonal weather patterns as well as in characteristics of hydroclimatic variables such as temperature and precipitation. Therefore, we selected CDs as the basis for our assessment of spatial hydrodynamic trends.
The USGS maintains a network of HCDN gages across the United States that represent streams with minimal or no anthropogenic disturbance or influence. As these streams are unimpaired by damming, artificial storage, and channel diversion for withdrawal and use, an analysis of their streamflow records allows assessment of hydrologic response to climate. As of 2009, there are 39 gaging stations in Texas that are continuously monitoring streamflow discharge (Lins 2012). Of these 39 stations, seven were selected for the present study based on the following criteria: (1) the gaging stations must be located on ecologically significant stream segments that serve a hydrologic function; and (2) daily streamflow records for the water year 1970 to water year 2019 (October 1, 1969 to September 30, 2019) must be available. The locations of these seven gages on their respective stream segments and the climate division they are contained within are shown in Figure 1. These seven gages span three CDs: site 0807000 is located in the Upper Coast (UC) CD, site 08171300 is located in the South Central (SC) CD, and sites 08165300, 08190000, 08190500, 08195000, and 08198500 are located in the EP CD. It must be noted that while site 08070000 is located in the UC CD, over 90% of its contributing watershed lies in the East Texas CD.  Daily discharge data were compiled for the seven gages using the EGRET (Exploration and Graphics for RivEr Trends) software package (Hirsch and De Cicco 2015) developed for use with RStudio (RStudio Team 2019). The FlowScreen package (Dierauer and Whitfield 2017) was used to develop and analyze baseflow statistics such as minimum, mean, and maximum baseflow at the gages for a user-defined period. The Eckhardt digital filter method (Eckhardt 2012) built into this package (which has been recommended by Xie et al. (2020) for the contiguous United States) estimates baseflow from streamflow discharge.
The watershed draining to each of the seven gages was first delineated. Then, the land use and land cover (LULC) characteristics and temporal changes of the seven watersheds were evaluated at 5-year intervals. The purpose of this exercise was to verify that the watersheds had undergone minimal LULC change, albeit over a 33-year timeframe, as only data from 1985-2017 were available from the USGS Land Change Mapping Assessment and Projection Datasets (USGS 2021). Nonetheless, validation of LULC changes further aids the attribution of hydrologic trends. If minimal or no LULC changes were present at the gages, we can attribute the hydrolog-ic trends observed to changes in climate (Lins 2012). Eight types of LULC are described in these datasets -developed, cropland, grass/shrub, tree cover, water, wetland, ice/snow, and barren. Ice and snow cover do not exist for any of the watersheds examined. The temporal trends in LULC for the seven watersheds are shown in Figure 2. It should be noted that only LULC for 5¬-year intervals beginning with the year 1985 are shown in this Figure. Developed land use is minimal (≤3%) and is only observed in watersheds draining to gages 08070000, 08171300, 08195000, and 08198500. In all four watersheds, there is no temporal change in the percent of the watershed area under developed use. The only watershed with any appreciable agriculture is that draining to gage 08070000. Cropland cover in this watershed shows very little change, ranging from 12% to 15% of the overall area depending on the time period of interest ( Figure 2a). Overall, with the exception of the watershed draining to gage 08171300 in the SC CD, where a slight increase (from 50% in 1985 to 60% in 2017) in grass/shrub cover occurred at the expense of tree cover, there was no notable change in LULC at any of the seven sites.
The location of the hydrologically unique stream segments ( Figure 1c) and the description of the functions they serve are shown in Table 1. All seven stream segments serve as sources of recharge for the respective aquifers they overlie. The East Fork of the San Jacinto River recharges groundwater to the Chicot Aquifer (which is part of the larger Gulf Coast Aquifer system), while all remaining segments overlie the Edwards Aquifer. Four gages (08171300 on the Blanco River, 08190000 on the Nueces River, 08190500 on the West Nueces River, and 08195000 on the Frio River) directly overlie the sensitive recharge zone of the Edwards Aquifer ( Figure 2). In addition to serving critical hydrologic roles, each of these seven stream segments is unique for one or more of the following reasons: (1) serving vital biological functions, (2) housing threatened or endangered species or unique communities, (3) providing high water quality or exceptional aquatic life use and aesthetic value, and (4) acting as a riparian conservation area (Table 1; TAC § 358.2 2012).
In addition to daily streamflow records, water level data from two groundwater wells (also referred to as "index wells") maintained by the EAA were compiled and included as part of the hydroclimatic assessment (EAA 2021b). As part of this monitoring effort, daily high water level data are available from two index wells: J17, representative of the "San Antonio Pool," and J27, representative of the "Uvalde Pool." Although daily-high data are available from 1932 for J17 and from 1942 for J27, only the daily-highs for the water years 1970-2019 were used in this analysis. Well J17 is located in the SC CD while J27 is in the EP CD ( Figure 3). Spring flows in the Aquifer help sustain seven endangered and one threatened aquatic species. Water withdrawals by pumping can detrimentally impact these flows and threatened species. Therefore, continuous monitoring of groundwater levels using these index wells is mandated. The EAA jointly uses the water level data from these wells and discharge data from two springs, the Comal Springs and the San Marcos Springs, to enforce groundwater withdrawal restrictions during periods of drought based on set criteria (EAA 2021b).
Monthly total precipitation, monthly total potential evapotranspiration (PET), and monthly average temperature were compiled from the University of East Anglia Climatic Research Unit (CRU)'s high-resolution gridded data, version 3.26 (Harris and Jones 2019). This dataset is presented at 0.5° x 0.5° resolution and has been widely used in catchment-scale studies (e.g., Demaria et al. 2013;Hajihoseini et al. 2015;Mahmood et al. 2019;Mutti et al. 2020). The weather data from observation stations reported by the National Climate Data Center (NCDC) were either discontinuous, sparse, or not available for part of our study area. As a result, the CRU dataset, which includes PET data, was used as an alternate source. Harris and Jones (2019) reported that while temperature and precipitation are primary variables based on observations, PET is a derived variable, computed from temperature, vapor pressure, and cloud cover. This dataset presents month-by-month variations in these climate variables for the period January 1901 to December 2017.
Data from the CRU grid that encompassed each watershed were compiled. In some instances, a watershed spanned two 0.5° x 0.5° grids; in these cases, climate data from the two grids were aggregated by area-weighted averaging. Data that were averaged in this fashion are still referred to in the singular (as "grid") for simplicity. The resulting pairing of gages and the CRU dataset is as follows: gage 08070000 is paired with the UC grid, gage 08171300 is paired with the SC grid, gage 08165300 is paired with the Edwards Plateau North (EPN) grid, and gage 08190500 is paired with the Edwards Plateau West (EPW) grid. The watersheds draining to gages 08190000, 08195000, and 08198500 are adjacent to each other and are all encompassed by the Edwards Plateau East (EPE) grid. The ncdf4 package (Pierce 2015) was used within RStudio to extract and analyze the precipitation, PET, and temperature datasets. Finally, Spearman's rank correlation was performed between streamflow depth (as calculated using Equation 1) and the three CRU climate variables to investigate the strength of their linear relationship. Spearman's "rho" or correlation coefficient was squared to give the coefficient of determination (referred to as r 2 in this study) to determine the variance in streamflow depth that can be explained by the climate variables. (1)

Hydroclimatic change indicators and metrics
To capture the magnitude and duration of the streamflow (low) extreme, the annual minimum of 1-day means and annual minimum of 30-day means were selected. These metrics represent the lowest single-day value in a water year and the lowest consecutive 30-day or monthly average occurring in that year. Considering that streamflow discharge measured at a gage comprises above-ground and below-ground components, averaging flows over monthly (or longer) time periods helps buffer short-term, or sudden, fluctuations and provides a means of analyzing the persistence of drier conditions. Additionally, a metric such as the 30-day minimum also represents a measure of a more prolonged hydrologic and, consequently, environmental stress. The next metric was the number of days below the low flow threshold, defined as the number of days in that year that the daily mean falls below the 25th percentile of the daily means of the entire period of the study (in this case, the water year 1970 to the water year 2019) at that location. The 25th percentile was adopted as the low flow threshold following the recommendations of The Nature Conservancy's Indicators of Hydrologic Alteration User Manual (2009). The aforementioned indicators are also a select subset of hydrologic alteration statistics prescribed by Richter et al. (1996).
In addition to these indicators, the runoff efficiency (RE) of the watershed was also computed as shown in Equation 2. RE is a measure of the fraction of precipitation that is converted to runoff and changes in this parameter may reflect climate variability (i.e., changes in temperature, precipitation, and PET). This metric was included in the current study as a tool for evaluating the nature of the precipitation-runoff relationship. To compute RE, the discharge data from the gages and precipitation data from the CRU dataset were paired in the same manner as described in Section 2.1. Lastly, the FlowScreen package was used to perform baseflow separation from the daily mean discharge data. The annual minimum, annual mean, and annual maximum baseflow were included as indicators in this study. (2)

Data Analysis
The modified Mann-Kendall test (MMK) proposed by Hamed and Rao (1998) accounts for autocorrelation by modifying the variance of the original Mann-Kendall test. It has been widely used in hydrologic studies for the detection of non-stationarity (e.g., Wahl et al. 2015;Venkataraman et al. 2016;Machiwal et al. 2019;Alashan 2020) and is employed in this study for the detection of monotonic trends. For the sake of brevity, the MMK test has not been discussed here (see Hamed and Rao 1998 for a comprehensive treatment).
The MMK was applied to the following hydroclimatic indicators: (1) for streamflow -the annual and seasonal 1-day minimum of means, the annual and seasonal 30-day minimum of means, the number of annual days below the low flow threshold, and the annual RE; (2) for baseflow -the annual minimum, the annual maximum, and the annual mean; (3) for groundwater levels -the annual minimum of, the annual mean of, and the annual maximum of daily-high water levels; and (4) for climate variables -the annual and seasonal mean temperature, the annual and seasonal total precipitation, and the annual and seasonal total PET. The significance of linear trends was assessed at p≤0.05. Additionally, the magnitude of the trends for the streamflow, baseflow, and groundwater levels were characterized using the Sen slope, an estimator which captures the linear rate of increase or decrease of a parameter over the time period of reference (Sen 1968). While the MMK helps ascertain whether a monotonic trend is present, the Sen slope helps compare the magnitude of this trend between different gages and climate divisions.

Drought Indices
Several standardized indices have been widely used in hydrological studies to investigate the length and severity of droughts. These include the Standardized Precipitation Index (SPI) (McKee et al. 1993), the Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al. 2010), the Standardized Runoff Index (SRI) (Shukla and Wood 2008), and others. The general procedure for the development of these indices involves identifying the probability distribution that best fits the time series aggregated over a period of interest (e.g., 1-month, 6-months, 12-months, etc.) and subsequently transforming this distribution to a normal distribution with zero mean and unit variance. In this study, the 12-month SPI (henceforth referred to as SPI-12) was computed for each of the five CRU grids. The SPI package in the R software environment (Neves 2013) was employed to compute the SPI-12. As the CRU dataset was available only until 2017, the SPI-12 was computed for water years 1970-2017. Values of this index that fall within -1 to +1 indicate "normal," or average, precipitation conditions, while those values that exceed +1 or are smaller than -1 indicate abnormally wet/above-average precipitation periods and abnormally dry/below-average precipitation periods, respectively. A similar procedure was followed for developing an index for streamflow, also referred to as the Standardized Streamflow Index (SSI) (following Vicente-Serrano et al. 2012), with zero mean and unit variance. The SSI was also computed on a 12-month scale and is hereafter referred to as SSI-12. For both the SPI-12 and the SSI-12, the number of months in each decade falling above-or below-average conditions was computed to compare drought severity.

Trends in climate variables
The trends in mean temperature, total precipitation, and total PET at both the annual and seasonal time scales for water years 1970 to 2017 were assessed using MMK (Table  2). Annual mean temperatures show significantly increasing trends in all five CRU grids. From a seasonal perspective, mean temperatures for all five grids for all four seasons are rising except the autumn and winter mean temperatures for the UC grid, as shown in Table 2. There are no significant trends in annual precipitation, but spring totals show declining trends in the EPE grid and the EPW grid. Finally, annual PETs show increasing trends in both the EPE and EPW grids. At these two grids, spring and summer PETs are also increasing. Additionally, autumn PETs show increasing trends in the EPN and EPE grids.
The strength of the linear relationship between annual streamflow depth in millimeters (mm) and the three climate variables, i.e., annual mean temperature, annual total precipitation (in mm), and annual total PET (in mm), is shown in Figure 4. For an explanation of the pairing of streamflow gages and CRU grids, please refer to Section 2.1; gage 08195000 was selected to represent the EPE grid. Annual streamflow is positively correlated with annual precipitation at all grids. The coefficient of determination for all five grids is statistically significant; the SC, EPE, and UC grids show the strongest cor-   relation between streamflow and precipitation, as indicated by the r 2 . As for PET, the annual streamflow in the SC, EPN, EPE, and EPW show a negative correlation, with the EPE showing the strongest r 2 of 0.52, indicating that just over half the variance in streamflow can be explained by PET. However, PET does not display a significant correlation with streamflow in the UC. Lastly, annual streamflow shows a negative correlation with annual mean temperature in the SC, EPN, EPE, and EPW, with the EPE again displaying the strongest r 2 of 0.37. Temperature seems to have no impact on the variance in streamflow in the UC, as shown in Figure 4.

Trends in streamflow and baseflow
The trends in annual and seasonal 1-day and 30-day minimum flows were assessed for the water years 1970-2019 at a significance level of 0.05 (Table 3). The UC CD (site 08070000) shows no significant trends in the 1-day or 30-day minimums at the annual or any of the seasonal scales. Likewise, the two westernmost sites in the EP CD (gages 08190500 and 08190000) show no trends at any time scale. It must be noted that gage 08190500 in the EPW frequently experiences several zero-flow days in a year. In 21 of the 50 water years included in this study, there were at least 60 days per water year with zero mean flow. However, analysis of daily mean flows at gage 08190000 did not reveal any zero-flow days for any of the water years chosen for the study. For other sites, the results are mixed. The SC CD (gage 08171300) shows significant declining trends in annual 1-day and 30-day minimum flows while, at the seasonal scale, summer 1-day and summer 30-day minimums show declining trends (Table 3). In the EP, sites 08165300, 08195000, and 08198500 all show declining annual 1-day and 30-day minimums. Site 08198500 shows significant declining trends for all four seasons of the year for all three low flow metrics while adjacent site 08195000 shows an identical pattern, albeit with autumn/fall (OND) 30-day minimums alone showing no significant trends (Table 3)

Low Flow Trends in Texas Stream Segments Serving Unique Hydrologic Functions
Texas Water Journal, Volume 14, Number 1 and OND 1-day and 30-day minimums, but there were no trends at this site in the winter season (JFM; Table 3). As for the magnitude of trends, despite showing significant declines in all annual and seasonal low flow metrics, site 08198500 had the lowest Sen slopes of all seven sites. In contrast, adjacent site 08195000 in the same EPE grid experienced the sharpest declines, with the OND and JFMs showing the largest slopes.
On an annual basis, the Sen slope for the 30-day minimum flows is larger than the 1-day counterpart at site 08171300, but on a seasonal basis, the Sen slopes for the summer 1-day and 30-day minimum flows are nearly identical. Trends in the number of days below the low flow threshold are shown in Figure 5. All sites except site 08070000 in the UC CD and site 08190500 in the EP CD show increasing trends ( Figure 5). The annual RE was computed for each of the five CRU climate grids, with gage 08195000 being representative of the EPE grid for water years 1970-2017. Patterns in annual RE are shown in Figure 6 and locations where this metric showed a statistically significant trend (using MMK at a significance level of 0.05) are indicated. At the EPW grid, which encompasses the westernmost site 08190500, the maximum RE observed was nearly 0.2, indicating that, at best, 20% of annual precipitation is translated to runoff (Figure 6e). The RE in the majority of the years is <0.1 in the EPW, which is explained by the number of low flow threshold days ( Figure  5). The largest interannual variability in RE is shown in the EPE grid ( Figure S2; violin plots showing the kernel density, median, and interquartile range are presented as supplementary material). Although the largest REs (slightly more than 0.6) were recorded at this site, this metric shows declining trends here as well as in the EPN grid (Figures 6c, 6d). Interannual variability in RE is also pronounced in the SC grid, albeit with no statistically significant trends (Figure 6a). The UC grid shows no significant trends either (Figure 6b). We note that the averages of the REs for the five grids over the chosen 48-year duration are similar to the long-term REs computed by McCabe and Wolock (2016) for the period 1951-2012 for the hydrologic units they fall within.
The MMK test for annual minimum, annual mean, and annual maximum baseflow showed no significant trends in these three metrics in the UC CD site (08070000) or in the two westernmost sites in the EP CD, i.e., gages 08190500 and 08190000 (Table 4). The three remaining sites in the EP CD (gages 08165300, 08190000, and 08190500), as well as the SC CD (gage 08171300), all show declining annual minimum baseflows. Two EP sites (08165300 and 08195000) show significant declining trends in annual mean and annual maximum baseflows (Table 4).

Trends in well water levels
Trends in the annual maximum, mean, and minimum of daily-high water levels recorded at wells J17 (SC CD) and J27 (EP CD) were assessed at a significance level of 0.05. There were no significant trends in any of the three water level metrics at well J17. However, the annual maximum, annual mean, and annual minimum of daily-high water levels all showed significant decreasing trends at well J27. On a comparative basis, the Sen slope of the annual minimum of daily-high water levels was larger than the annual mean and annual maximum of daily-high water levels at this well.  From Figure 3, it is evident that well J27 lies in the artesian zone downgradient of gages 08190000 and 08195000, as well as in the same EPE grid as these two gages and gage 08198500. Water levels here are possibly influenced by the hydrology of the upgradient streams (losing streams) located in the recharge zone as well as climate variables such as temperature, precipitation, and PET. To further explore these relationships, cross-correlations and Spearman correlations were performed to estimate the r 2 . The cross-correlation functions (CCFs) between the annual mean of the daily-high water level (AMDHWL) at well J27 and the three climate variables (on an annual scale) are shown in Figures 7a-7c. While precipitation appears to have no cross-correlation with the AMDHWL (Figure 7b), both temperature and PET show negative CCFs with a lag of approximately 1-2 years, indicating that above-average annual temperatures and PETs precede below-average groundwater levels by approximately 1-2 years. An even stronger CCF is found between annual mean streamflow at the upgradient gages (08190000 and 08195000) and the AMDHWL, as seen in Figures 7d and 7e. The CCFs are positive and peak at a 1-year lag, suggesting the influence of recharge to well J27 from above-ground flows at these two gages. A similar pattern is evident with the annual mean baseflows at these two gages (see Figures 7f and 7g). The strength of the linear relationship between the AMDHWL and the related hydrological variables is demonstrated using the r 2 metric in Figure 8. The negative correlation between the AMDHWL and both annual mean temperature and annual total PET are evident from Figure 8; the r 2 for both climate parameters is roughly the same (0.18) and is statistically significant. Weak positive correlations between the AMDHWL and annual mean streamflows as well as baseflows are also evident (Figure 8). The strongest r 2 occurs at gage 08195000 upgradient of well J27; 32% and 40% of the variance in AMDHWL are explained by streamflow and baseflow, respectively.

Analysis of drought indices
The percent of each decade spent under abnormally dry (i.e., below-average) and abnormally wet (i.e., above-average) conditions according to the SPI-12 and 12-month Standardized Streamflow Index (SSI-12) was computed and is shown in Figure 9. Abnormal conditions are defined as periods when the index exceeds +1 or drops below -1, while normal, or "average," conditions are characterized by values of the index between -1 and +1. Figure 9b shows that very little of each decade leading up to the 2010s was characterized by below-average flows in the EPN and EPE grids. However, more than 66% of the 2010s (water years 2010-2017) were characterized by below-average flows at these two grids. This observation is in stark contrast to the trends in the SPI-12 for these two grids (Figure 9a), where only 15% or less of the 2010s are classified as drier-than-average, suggesting a pronounced impact of temperature-influ-enced drying. Another interesting observation is the lack of above-average flow periods at the gages in the EPN and EPE during the 2010s, despite experiencing wetter-than-average conditions at least 15% of the time (Figure 9b and 9a, respectively). The temporal variations in the SPI-12 and the SSI-12 are shown as supplementary material ( Figure S4).

Hydroclimatic trends, drivers and implications
At site 08070000 located in the UC CD, it appears that the increasing trend in annual mean temperatures has not impacted PET or streamflow. Although this gage is located in the UC CD, its contributing drainage basin is almost entirely located in the adjacent East Texas CD. Jiang and Yang (2012), Venkataraman et al. (2016), and Crawford et al. (2019) note that the eastern extreme of Texas (which encompasses the watershed contributing to gage 08070000) benefits from moisture surplus due to proximity to the Gulf of Mexico and generally suffers milder droughts relative to the rest of the State.
In the SC CD, statistically significant increases in annual and seasonal mean temperatures were detected at the SC grid, yet there were no trends in precipitation or PET. This gage (08171300) overlies the Edwards Aquifer recharge zone, where low flows and baseflows have been decreasing over the 50-year period beginning in 1970. Nearly one-third of the 2010s was spent in below-average SSI conditions -the worst among the five . However, in comparison with other CDs, this site endured the fewest below-average flow months in the 2010s. As such, the decrease in low flows and RE at this location may be temperature-driven, but modeling studies that explicitly account for the influence of temperature on streamflow are needed for further validation.
For the EP CD, generalizations cannot be made about any hydroclimatic variables except that annual and seasonal mean temperatures show increasing trends. While annual 1-day minimum and 30-day minimum flows exhibit decreasing trends at gage 08165300 on the South Fork of the Guadalupe River (in the EPN) and gages 08195000 on the Frio River and 08198500 on the Sabinal River (both in the EPE), gage 08190500 on the West Nueces River (EPW) shows no trends whatsoever; the same patterns are exhibited in baseflow. It appears that the West Nueces River in the EPW is intermittent, experiencing many zero flow days in a year. Consequently, there are no trends in any of the streamflow metrics or baseflow (which was separated from the streamflow data) despite increasing annual mean temperatures and PET. In fact, long periods of zero flow days, some lasting five consecutive months, have been reported in the West Nueces River (Thomas et al. 2019;Hackett 2019).
The drying trend in the EPN and EPE is further evident in the increasing number of days below the low flow threshold and decreasing runoff efficiencies. The EPN and EPE grids experienced no wetter-than-average months in the 2010s, and of the five chosen grids, experienced the worst streamflow droughts, as indicated by the percentage of the 2010s in below-average SSI conditions. The drying pattern in the EPE is particularly worrisome considering its moderate influence on daily-high water levels in well J27 located in the same climate grid. Lindgren et al. (2004) reported no long-term declines in water levels in the Edwards Aquifer, but their study was limited to the 20th century. They mentioned that water levels showed rapid recoveries after periods of drought and that the highest water levels were observed in the 1990s. Although we included only two wells in the Edwards Aquifer as part of our study, one of which (J17) showed no trends in daily-high water levels, there is evidence that daily-highs in well J27 have been decreasing since 1969. From a natural resource management perspective, water levels in the J17 and J27 wells are used by the EAA as the criterion for distinguishing stages of drought as part of a critical period management (CPM) plan. In Uvalde County, the water level of well J27 has been reported to be the most suitable indicator of drought severity by Green and Bertetti (2010), as river discharge did not appear to be useful. However, the cross-correlation at a 1-year lag between mean flows at gages 08190000 and 08195000, and the mean of daily-high water levels at well J27, as well as the moderately strong but statistically significant r 2 between the two, suggest that discharge at these two gages merits consideration for an early warning or preliminary drought trigger system for Uvalde County.
Although annual minimum flows and annual baseflows are decreasing in the EPN and EPE, the magnitude of these declines is higher in the EPE. This difference between the two grids is likely due to the combined effect of increasing temperature and increasing PET in the EPE as opposed to increasing temperature alone in the EPN. It is also worth noting that spring precipitation is decreasing in the EPE. Precipitation regimes in the EP CD are generally bimodal, with spring (May) and end-of-summer (September) accounting for the majority of annual precipitation. Therefore, decreasing spring precipitation possibly has a role in the drying observed here. Thomas et al. (2019) found similar drying patterns in streamflow in parts of the Upper Nueces River Basin, portions of which overlap our study area. They did not find any conclusive trends in ET in what is essentially the EPE grid of our study area for the period 1970-2015. However, their findings are based on ET data developed on a continental scale. Using the CRU climate dataset, which allows for analysis at a more localized scale, we have found increasing PET trends in the EP. The results of the correlation analyses also showed that the EPE grid had the highest r 2 between streamflow and PET as well as temperature of the five chosen climate grids. These r 2 (0.52 for PET and 0.37 for temperature) suggest the moderate influence of these meteorological parameters on streamflow. Considering the minimal anthropogenic impact on these watersheds, the stream drying observed here may be climate-driven.
As for secondary factors involved in the drying pattern observed in the EPN and EPE, we first consider the karst landscape of the EP and its impact on rainfall-runoff relationships. Wilcox et al. (2008) reported that where soils are shallow and underlain by impermeable limestone in this region, overland flow dominates subsurface flow. Additionally, they highlighted the presence of overland flow zones in areas with certain types of vegetation, i.e., woody plants versus grass and shrub cover. Wilcox and Huang (2010) further suggested that degradation of karst landscapes may result in declines in groundwater recharge and, subsequently, baseflow, but above-ground river flows may recover with an increase in woody plant cover. The implications of these two studies are that in the EP, particularly the EPE grid, spatial variations in the karst landscape may result in (a) a greater fraction of overland flow versus subsurface flow, which may lead to greater exposure to the elements in a drying climate; or (b) reductions in baseflow where these landscapes may be degraded. Lindgren et al. (2004) emphasized that the recharge zone of the Edwards Aquifer is characterized by a "dual-or triple-porosity nature" and that groundwater flow in this zone is poorly understood. Recent studies, such as Kromann (2015) and Hackett (2019), report that the nature of streamflow losses and gains in these segments and their subsequent role in drying patterns are unclear and merit further research. The last factor that cannot be ignored is the role of surface water governance. Although groundwater use in the Edwards Aquifer is strictly regulated by the EAA, all surface water in this watershed is owned by the State of Texas and appropriated to users through a system of water rights permits. Thomas et al. (2019) noted that withdrawal of surface water beyond the allocated quota may occur in parts of the EP CD (which fall under the EAA's jurisdiction), particularly in times of low streamflow. Such violation of the honor system (which involves self-reporting of water extraction) may go unnoticed and unreported. It must be added that the Blanco River in the SC climate division may suffer from the same exploitation of water rights since it does not fall under the purview of a Watermaster system.

Ecological implications
Regarding ecological implications, the reduction in baseflow we note in our study may influence stream communities of the EP and may provide a preliminary indication of a changing flow regime. Reduced baseflow and drying events/droughts have demonstrated effects on the communities of other streams. For example, even though the low presence of riffles in streams with intermittent flows can lead to lower species richness than in perennial streams, macroinvertebrate communities are still found to be diverse (Santos and Stevenson 2011). Inter-estingly, the communities of streams with intermittent flows can be distinct from perennial streams in terms of functional feeding groups and benthic communities (Santos and Stevenson 2011), rare and endemic niche specialists , and a combined biodiversity composed of aquatic and terrestrial species due to the dry-wet regime of the system (Bunting et al. 2021). Streams with intermittent flows have diverse and shifting communities during lentic (flowing), lotic (ponding), and dry phases such that they contribute greatly to the overall biodiversity of the entire catchment, both aquatic and terrestrial Hill and Milner 2018). In the EP, hydrological connectivity of the streams is critical to the resilience of the basin community, as baseflow is apt to decline as drought worsens. Moreover, the dry-wet phases of a stream require that both aquatic and terrestrial communities be characterized and considered in management, as they are both affected and can colonize a streambed quickly in either flow regime (Bogan et al. 2017;Bunting et al. 2021). This suggests that even terrestrial species of ecological concern in EP, such as the Texas Snowbell (Styrax platanifolious ssp. texanus), could potentially be influenced by changing flow regimes.
The management and conservation of ecologically significant stream segments ideally focus on maintaining ecological resilience (here defined as the ability of an ecosystem to maintain structure and function in the face of disturbance, per Holling 1973). The findings of this study agree with those of others that climate exerts a high-order environmental control on low-and no-flow stream conditions. However, regional-scale climatic, physiographic, and anthropogenic factors play important roles in determining the flow regime of streams (Reynolds et al. 2015;Hammond et al. 2020). Drought protection of ecologically significant stream segments should consider habitat diversity to preserve ecological functions. For example, perennial pools and flowing reaches provide drought refuges and habitat for newly colonizing taxa Hill and Milner 2018), and headwaters, especially in forested catchments, host high biodiversity (Storey et al. 2011;Bogan et al. 2015). In general, knowledge of the spatial distribution of perennial and intermittent river channels in a river basin would optimize such protection plans (González-Ferreras and Barquín 2017). Furthermore, an ecological understanding of the life history traits of the organisms relying on drought refuges, especially endangered species , would help provide targeted conservation management plans for ecologically significant stream reaches.

CONCLUSION
Investigation of low flows is a critical part of evaluating the overall health of a stream and is imperative for long-term natural resource management. In this study, we have used a vari-ety of metrics to assess the low flow characteristics of stream segments in Texas that serve a unique hydrologic function for the period covering water years 1970 through 2019 using discharge data from seven USGS gages. Although annual mean temperatures have been increasing in all climate divisions chosen in this study, critical inter-and intra-climate division differences highlight the spatially diverse nature of the State of Texas. As such, there are no significant streamflow trends in the gage located in the UC CD, the watershed for which lies in the East Texas CD, likely due to proximity to the Gulf of Mexico. The Blanco River in the SC CD, and the trio of the South Fork of the Guadalupe River, the Frio River and the Sabinal River in the EP CD have undergone pronounced drying since the water year 1970. We performed LULC analysis of the watersheds contributing to these gages and confirmed that they have undergone little to no change over time, indicating minimal anthropogenic influence. Therefore, the results of the correlation analysis with climate variables and the comparison of drought indices suggest that the drying we have observed may be climate-driven. It must be noted that our findings and interpretations are based on a small subset of stream gages confined to three climate divisions in the State, and therefore far-reaching conclusions or generalizations about regional patterns or trends cannot be made.
Overall, in planning for changes associated with climate change and human water demand, several studies concur that intermittent streams are critical components of a river basin in terms of biodiversity, community dynamics, biogeochemical cycling, ecosystem services, and ecological resilience, even during dry phases. This study finds evidence of increasing temperature, declining spring precipitation, increasing PET, and declining minimum flows and daily-water level highs for the EP, although there are seasonal and physiographic differences in these trends among sites. In general, if the EP were to experience more temperature-driven drying and drought conditions in the future, the streams would provide important drought refuges at perennial pools and hyporheic zones in addition to the habitat and ecosystem services they already provide as flowing waterbodies. Moreover, the role of streams in recharging the Edwards Aquifer will remain important even when the streams do not have surface flow. The EP streams investigated in this study will likely remain ecologically significant in the future and the use of climatic and land use variables to monitor and predict conditions in the region will be critical for their conservation and management. LULC codes: 1 -developed land; 2 -cropland; 3 -grass/shrub cover; 4 -tree cover; 5 -water; 6 -wetland; 8 -barren; double digit codes show change from one category to another. For example, 12 indicates change from class 1 (developed) to 2 (cropland) during that 1-year time frame.    Table S1. Summary of trends and changepoints in annual minimum, mean, and maximum of daily-high water levels at wells J17 and J27 (significance assessed at p≤0.05). Sen slope values are shown in parentheses.