*Article* **Temporal E**ff**ects of Groundwater on Physical and Biotic Components of a Karst Stream**

**Tao Tang 1, Shuhan Guo 1,2, Lu Tan 1, Tao Li 2,3, Ryan M. Burrows <sup>4</sup> and Qinghua Cai 1,\***


Received: 11 May 2019; Accepted: 19 June 2019; Published: 21 June 2019

**Abstract:** Although most lotic ecosystems are groundwater dependent, our knowledge on the relatively long-term ecological effects of groundwater discharge on downstream reaches remains limited. We surveyed four connected reaches of a Chinese karst stream network for 72 consecutive months, with one reach, named Hong Shi Zi (HSZ), evidently affected by groundwater. We tested whether, compared with other reaches, HSZ had (1) milder water temperature and flow regimes, and (2) weaker influences of water temperature and flow on benthic algal biomass represented by chlorophyll *a* (Chl. *a*) concentrations. We found that the maximum monthly mean water temperature in HSZ was 0.6 ◦C lower than of the adjacent upstream reach, and the minimum monthly mean water temperature was 1.0 ◦C higher than of the adjacent downstream reach. HSZ had the smallest coefficient of variation (CV) for water temperature but the largest CV for discharge. Water temperature and discharge displayed a significant 12-month periodicity in all reaches not directly groundwater influenced. Only water temperature displayed such periodicity in HSZ. Water temperature was an important predictor of temporal variation in Chl. *a* in all reaches, but its influence was weakest in HSZ. Our findings demonstrate that longer survey data can provide insight into groundwater–surface water interactions.

**Keywords:** groundwater–surface water interactions; time series; water temperature regime; flow regime; chlorophyll *a*; wavelet analysis; convergent cross-mapping

#### **1. Introduction**

Most streams and rivers are groundwater-influenced systems [1,2]. Groundwater commonly upwells into surface waters at the hyporheic zone [3,4], which influences the physical and biotic characteristics of surface waters substantially [5–8]. The importance of groundwater to surface lotic ecosystems has been documented extensively, with most of the work focusing on distinct functions of the upwelling zone. Upwelling zones have been observed to have more favorable conditions for many organisms, including stable temperature and flow conditions as well as higher concentrations of dissolved nutrients [9,10]. Consequently, upwelling zones generally harbor higher biodiversity and more specialized biotic assemblages compared with other lotic environments [11,12]. The upwelling zone also serves as a refuge for lotic organisms to survive adverse environmental conditions such as extreme temperatures, floods, and stream drying [13,14]. A handful of studies have estimated the ecological impact of groundwater discharge for organisms in downstream reaches; however, findings are equivocal. Some studies demonstrate that important ecological characteristics, such as

primary production and biodiversity, of downstream reaches benefit from groundwater inputs [15]. Further, coldwater fish have been confirmed to use groundwater-fed reaches as a thermal refuge during warming periods [16,17]. In contrast, other studies did not find significant groundwater effects for aquatic organisms [18,19].

Study duration is likely a major reason underlying variation in the reported reach-scale physical and biotic effects attributed to groundwater discharge. Most studies have evaluated groundwater effects over short time scales (usually never exceeding two years) and without the adequate characterization of the changes in seasonal patterns and processes [2,16,18,20]. Since the influence of groundwater for lotic ecosystems may vary remarkably among seasons and years [21,22], inconsistent results may be expected if study designs do not adequately account for this temporal variation. Moreover, short-term studies often cannot establish accurate or meaningful quantitative relationships between measured variables due to limited data. This limitation hampers mechanical exploration of pathways describing how groundwater influences surface waters. Therefore, studies of longer duration and with a more intensive survey effort are needed for exploring long-term and quantitative effects of groundwater discharge on downstream ecosystems, which will be valuable for a better understanding of groundwater–surface water interactions.

Short-term studies may also fail to accurately characterize important groundwater effects on stream physicochemical conditions. For instance, stable water temperature and flow conditions are key features of groundwater-influenced stream reaches, and these conditions may mitigate temporal fluctuations of water temperature and flow in downstream reaches [23]. However, due to short durations, many studies only measure changes in the magnitude of downstream water temperature or flow (i.e., the incremental deviation from the average conditions caused by groundwater inputs) [24,25]. Importantly, other widely used regime descriptors for any continuous phenomena, including variability (temporal fluctuations across regular time intervals), frequency (the number of events occurring at a certain magnitude or exceedance), timing (dates that particular events occur), and predictability (the reliability of event recurrence), have been poorly documented [26,27]. Consequently, it is unclear if or how groundwater can change other aspects of temperature and flow regimes other than the response magnitude in downstream reaches. Such knowledge has ecological value because the magnitude, variability, frequency, timing, and predictability are all important aspects influencing lotic ecological structure and function [28,29].

In the present study, four connected reaches in a Chinese Karst stream were surveyed continuously for 72 months, one of which is known to receive direct groundwater inputs. We monitored water temperature, discharge, and benthic algal biomass (represented by Chlorophyll *a* [Chl. *a*] concentrations) to evaluate how groundwater discharge affects these important physical and biotic components of lotic ecosystems. Benthic algae are a high-quality food source in streams and can be the main energy source to the aquatic food web [30,31]. Water temperature and discharge are commonly influenced by groundwater discharge [9,10] and are also major drivers of the temporal dynamics of benthic algal biomass [32,33]. Our main objective was to characterize the water temperature and discharge regimes for each reach. Further, we assessed reach-specific relationships between benthic Chl. *a* concentrations and water temperature as well as discharge. We tested two hypotheses: the reach most affected by groundwater will have (1) milder water temperature and flow regimes and (2) weaker influences of water temperature and flow on benthic algal biomass dynamics. Support for these two hypotheses will provide evidence that groundwater has an important influence on temporal dynamics of adjacent downstream reaches in the karst stream.

#### **2. Materials and Methods**

#### *2.1. Study Region*

We conducted the study on the Nan River, located in Shengnongjia (31,015 –31,075 N, 109,056 –110,058 E), the western mountainous region in Hubei province of China. Elevation in this region varies from 420 to 3015 m (with a mean elevation >1500 m), resulting in an altitudinal

climate gradient. The annual mean air temperature ranges from 7.1 ◦C in the western region to 14.5 ◦C in the eastern region, with the coldest and warmest months in January and July, respectively [34]. The annual precipitation commonly varies between 800 and 2500 mm, in which more than 85% of rainfall occurs in April to October with December to February accounting for <8% [34]. Karst terrain is a typical geomorphological characteristic in Shengnongjia where ~30% of bedrock is limestone. Due to the high solubility of this type of rock, there are enormous underground rivers, karst caves, swallow holes, groove, and springs in this region [34,35].

We surveyed four reaches in the headwaters of Nan River (Figure 1). The streambed substrate is dominated by gravels, pebbles, and cobbles, with bedrock of limestone [34]. Sheng Nong Yuan (SNY) and Da Long Tan (DLT) are two first-order reaches with elevation of 2320 m and 2220 m, respectively (Table 1). SNY is located in the mainstream, while DLT is in an adjacent tributary (Figure 1). Swallow holes are common along the 1-km-long stream channel just downstream of the confluence of DLT tributary (Figure 1). Surface water seeps into swallow holes, leading to reduced flow in the reaches immediately downstream. The downwelling reaches (i.e., sections with swallow holes), which were not surveyed in this study, have been observed to dry up from late December to February when rainfall is minimal (personal observation). Groundwater recharges to the surface channel through discrete springs in reaches downstream of the swallow holes and ensure sufficient base flow to maintain surface water during dry periods. Hong Shi Zi (HSZ) is a second-order downstream reach approximately 2 km downstream from the confluence point of DLT with an elevation of 1950 m. Several visible springs can be observed along this reach; indicating that this reach is influenced by groundwater discharge. Ji Zi Gou (JZG) is another second-order reach about 4 km downstream HSZ, located at an elevation of 1820 m. All the four study reaches have permanent surface flow. The DLT tributary reach had the narrowest wetted width and the shallowest water depth during the study (Table 1). As for the three mainstream reaches, mean wetted width and water depth both increased longitudinally among reaches from SNY to JZG. All the four reaches had closed riparian canopy and low degree of human disturbance because of effective protection of this region.

**Figure 1.** Locations of the 4 stream reaches in the headwaters of Nan River. The ellipse circles illustrate stream channel where discrete swallow holes and springs were observed. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.


**Table 1.** Channel morphology characteristics for the four study reaches.

#### *2.2. Data Collection and Preparation*

The four reaches were surveyed monthly from July 2011 to June 2017; however, data was not collected in April and July of 2015 due to large floods. We monitored stream discharge, air temperature, and water temperature for each reach on all sampling occasions. Mean reach discharge was calculated from measurements at three randomly selected cross-sections along a 50-m reach. At every 50-cm interval along each cross-section, flow velocity was measured using a Global Water Flow Probe Hand-Held Flowmeter (Xylem Inc., White Plains, NY, USA) and water depth was recorded using a tape. The wetted width of each cross section was also measured with a tape. Discharge (m3 s<sup>−</sup>1) was calculated using the velocity-area method for each cross-section [36]. Air and water temperature were measured and averaged from the same three cross-sections with an YSI Pro Plus multimeter (YSI Inc., Yellow Springs, OH, USA). Although we did not sample all the reaches at the same time on each occasion, fieldwork in all reaches was completed within a 2 to 3-hour period to minimize variation due to time of sampling.

Benthic algae were sampled at the same aforementioned three cross-sections for each reach during each occasion. Five stones of 10–20 cm diameter were randomly selected at each cross-section, with a total of 15 stones for each reach. We defined the sampling area for each stone with a circular lid of 2.7 cm radius. The stone surface within the lid was vigorously scrubbed using a nylon brush and rinsed 3–4 times with distilled water. All subsamples from the same cross-section were combined into one composited sample and the volume was recorded [37]. Algal sample was passed through a glass fiber filter (Whatman GF/F) in the field and kept dark and cool. In the laboratory chlorophyll *a* was extracted with 90% acetone and absorbance values at 630, 645, 665 and 750 nm were measured using a Shimadzu UV-1601 spectrophotometer (Shimadzu Corp., Kyto, Japan), following a standard analysis procedure [38]. Benthic Chl. *a* concentration for each cross-section was calculated as mg per stone surface area (mg m−2). Chl. *a* concentration for each reach was averaged from the three cross-section replicates.

#### *2.3. Statistical Analysis*

Although our data did not allow us to directly quantify variation in groundwater-surface water exchange among the study reaches and over time, we made use of long-term (six years) measurements to infer spatial and temporal variation in the relative influence of groundwater. This approach has been employed elsewhere [7,16,19,25]. We then relate this understanding of the relative variation in groundwater influence to the observed spatial and temporal variation in the physical and biotic parameters that we assessed.

We performed bivariate linear regression by using water temperature as the response variable and air temperature the predictor to investigate the degree of groundwater influence for each reach. Since monthly water temperature data is not generally autocorrelated, simple linear regression is effective in modeling such a relationship [28]. We expect that the reach most influenced by groundwater (i.e., HSZ) would have a greater intercept and a lesser slope for the regression equation than reaches not influenced by groundwater [39].

We calculated nine metrics characterizing stream water temperature and discharge regimes during the study period for each reach. These metrics represent the water temperature and discharge regime magnitude, variability, frequency, and timing. Magnitude was described with the total mean value, the maximum monthly mean value (Mean\_max), and the minimum monthly mean value (Mean\_min) of water temperature and discharge. Variability was represented by the range and coefficient of variation (CV) of the two variables. We counted the number of months with water temperature >10 ◦C for delineating temperature frequency because 10 ◦C is an important temperature threshold for many aquatic metabolic processes [40,41]. For instance, algal biomass was greatest when water temperature was at or below 10 ◦C in several forested streams [41]. We also selected 3 ◦C as another threshold for calculating water temperature frequency because we observed that the study reaches displayed substantial differences among the reaches when temperature was lower than this value. For discharge frequency, we counted number of months with discharge <0.02 or >0.20 m<sup>3</sup> s−<sup>1</sup> which represents approximately the 25th and 75th percentiles of the distribution of discharge across all reaches. For the timing metrics, we identified specific months when the maximum (M\_mean\_max) and minimum (M\_mean\_min) monthly mean water temperature or discharge was reached. We also calculated all the aforementioned nine metrics for Chl. *a*, because algal biomass also varies temporally [42].

We used wavelet analysis to describe regime periodicities in water temperature, discharge, and Chl. *a* for each reach. Wavelet analysis is especially effective in detecting and differentiating scale-specific dynamics from noisy, complex and usually nonstationary time series data that violates assumptions of many other statistical methods [43,44]. This method performs a local time-scale decomposition (i.e., local wavelet transform) of the time series by estimating its spectral characteristics as a function of time [43]. Significance of local wavelet spectrum was tested by comparing against a background red noise spectrum (representing low-frequency/long-period cycles of the time series with autocorrelated residuals) simulated with Monte Carlo procedure [45,46]. Wavelet analysis was conducted by using 'wt' function in R package 'biwavelet', with continuous Morlet wavelet transforming as the base function due to its good balance in time and frequency localization [45]. To further determine which periodicities were most important throughout the study duration, we calculated time-averaged wavelet spectrum for temperature, discharge, and Chl. *a* with 'analyze.wavelet' function in package 'WaveletComp'. For each site, there was a total of 3–4 missing data during the study period. Since missing data are not allowed in wavelet analysis, for each time series, we interpolated missing values with the mean of all available data from the same month adjusted for that year's mean value using the decomposition method described by Cloern and Jassby [42,47].

We conducted convergent cross-mapping (CCM) to identify influences of water temperature and discharge on temporal fluctuations of Chl. *a*. This novel nonparametric method has been described in detail elsewhere [48–50], so we provide only a brief introduction here. Based on state space reconstruction, CCM measures the skill of cross-mapping (i.e., prediction) between two variables (e.g., x and y) as test for causality. Following Takens' Theorem, it is assumed that if x influences y, then y would contain unique information from x. Therefore, historical values of x can be reconstructed from y alone. Conversely, x does not contain information from y, hence will be incapable of cross-mapping y [50]. To accomplish the test, Simplex projection with a time delay embedding from y is first used to predict historical values of x [48]. Then, Pearson's correlation coefficient (ρ) between predicted x and measured x values is calculated to indicate the ability of y cross-mapping x. The ability of x cross-mapping y is also estimated with similar procedures. It is determined that x influences y if ρ value for y cross-mapping x is higher than that of x cross-mapping y, and vice versa [51]. Significance of the result was tested with Ebisuzaki phase-shift null model. This randomizing procedure removes dynamic correlations between time series, but preserves most of short-term behaviors [52]. A CCM result is considered significant if the ρ value for x cross-mapping y exceeds the range of the 5th and 95th percentiles of corresponding estimates for the null model. In the present study, we estimated reach-specific effects of stream water temperature and discharge on Chl. *a*. If temporal dynamics of Chl. *a* was driven by water temperature or discharge, Chl. *a* would be more effective in cross-mapping

the other two variables than using the other two variables to cross map Chl. *a*. We applied simplex projection to calculate the optimal embedding dimension for individual time series, and set library size varying from 1 to 60 in steps of 6 with 100 bootstrapped library samples. CCM was performed with time lags of ±3 months to account for possible time delayed effects [53]. For the time lag with the highest cross map skill, we further tested the significance of the results by generating a null distribution with 100 surrogate data. All the time series were first-differenced to remove any trends, and normalized to zero mean and unit variance for allowing time series comparisons [48]. CCM was performed with package 'rEDM'.

Additionally, to investigate whether contrasting seasonal hydrological conditions may influence water temperature, discharge, and Chl. *a,* we calculated the mean values (± standard deviation) of these parameters for each reach during a typically low-flow (December to February) and high-flow (May to July) period in this region [34]. We assessed whether these parameters differed among reaches within seasons using Kruskal–Wallis test by ranks, with pairwise comparisons by Mann–Whitney tests with Bonferroni correction.

#### **3. Results**

#### *3.1. E*ff*ects of Groundwater on Stream Water Temperature and Discharge Regimes*

Reach-specific regression relationships between air and stream water temperature were detected (Table 2). HSZ had the highest intercept of 4.162 and the shallowest slope of 0.352 for the relationship between air and stream water temperature. HSZ also had a distinct water temperature regime when compared with other reaches (Table 3). Mean water temperature increased from SNY to JZG as elevation decreased. However, the maximum monthly mean water temperature (Mean\_max) in HSZ was 0.6 ◦C lower than the adjacent higher elevation reach of DLT (Table 3). Further, the minimum monthly mean water temperature (Mean\_min) was 1.0 ◦C higher than the adjacent lower elevation reach of JZG (Table 3). Moreover, HSZ had a narrower temperature range than that in DLT and had the smallest coefficient of variation (CV) among all the four reaches (Table 3). In accordance with low temperature variability, there was only two months in HSZ when water temperature was <3 ◦C (Table 3). In comparison, for JZG, where mean water temperature was 1.5 ◦C higher compared to HSZ, mean monthly water temperature was <3 ◦C on eight occasions (Table 3). The frequency for months with water temperature >10 ◦C in HSZ was similar to that recorded in DLT (Table 3). Although influenced by groundwater, HSZ had same timings for the maximum and minimum monthly mean values of water temperature as that in other reaches. All reaches also had the same temperature periodicity that ranged from 10 to 14 months, with the 12-month periodicity having the most significant wavelet spectrum power across the study duration (Table 4, Figure 2). That is, stream water temperature had an annual cycle that peaked in August and was lowest in January (Figure 2). During the high flow period, the mean water temperature in HSZ was similar to that in the higher elevation reach of DLT (~10 ◦C for both reaches). Moreover, HSZ had a higher mean water temperature than that in the lower elevation reach of JZG (3.9 vs. 3.5 ◦C) during the low flow period, despite the difference was not significant (Appendix A Figure A1).

The discharge regime in HSZ was also different from that recorded in the other three reaches. The discharge magnitude and variation range generally increased with stream size (Table 3). The mean and range of discharge in the largest reach (JZG) was approximately 6 and 5.5 times greater, respectively, than that recorded in the smallest reach (DLT) (Table 3). However, other regime metrics did not display such stream size-related trends. HSZ had the highest discharge CV, and the highest frequency of months with discharge lower than 0.02 m<sup>3</sup> s−<sup>1</sup> (Table 3). Discharge displayed a significant 12-month periodicity in the all reaches across the study duration with the peak in May and trough in January or February (Table 4, Figure 3). However, significance of such periodicity (indicated by black line circles) varied among years and among reaches (Figure 3). Besides the annual periodicity, discharge also displayed a significant 2–6 months periodicity in 2013 across all reaches (Figure 3). During the low flow period, mean discharge in HSZ was 0.03 m<sup>3</sup> s<sup>−</sup>1, which was lower than the sum of mean discharge of the two upstream reaches (SNY: 0.03 m<sup>3</sup> s−1, DLT: 0.02 m3 s−1). In contrast, this pattern was not evident during the high flow period (Appendix A Figure A1).

**Table 2.** The Y intercept, slope, and adjusted R<sup>2</sup> (R2 *adj*) for the regression between stream water temperature and air temperature from July 2011 to June 2017 for each study reach. All regression models are significant with *p* < 0.001. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.


**Table 3.** Values of stream water temperature, discharge, and benthic chlorophyll *a* concentration (Chl. *a*) regime metrics for each study reach during the study. Mean\_max = the maximum monthly mean value; Mean\_min = the minimum monthly mean value; M\_mean\_max = the month with the maximum monthly mean value; M\_mean\_min = the month with the minimum monthly mean value. CV: coefficient of variation. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.


**Table 4.** Significant (*p* < 0.05) periodicities (Unit: month) in stream water temperature (WT), discharge, and benthic chlorophyll *a* concentration (Chl. *a*) across July 2011 to June 2017 identified by comparing time-averaged wavelet spectrum to a null red-noise power spectrum by wavelet analysis for each reach. Numbers in parentheses represent the periods having the highest average wavelet power. -: no significant period is detected. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.


**Figure 2.** Continuous wavelet power spectra showing the strength of the periodicities of stream water temperature changes between July 2011 and June 2017 for each stream reach. The colors code for power intensity from dark blue (low intensity) to dark red (high intensity), and solid black line circles indicate significant (*p* < 0.05) coherent time–frequency regions. The shaded area is the time–frequency region affected by edge-effects and was not included in the analyses. The original time series are demonstrated above each wavelet plot. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

**Figure 3.** Continuous wavelet power spectra showing the periodicity of discharge from July 2011 to June 2017 for each stream reach. Further information on wavelet spectra are described in Figure 2. The original time series are demonstrated above each wavelet plot, with different value ranges to illustrate reach-specific variation. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

#### *3.2. E*ff*ects of Groundwater on Associations between Chl. a Dynamics, Water Temperature, and Discharge*

Across all stream reaches, the range of Chl. *a* concentrations and frequency of months with Chl. *a* concentration higher than 10 mg m−<sup>2</sup> increased as elevation decreased (Table 3). The temporal dynamics of Chl. *a* concentrations in HSZ was, however, largely different than in the reaches less affected by groundwater (Table 3). Specifically, Chl. *a* magnitude was lower in HSZ when compared with the smaller upstream reach, SNY (Table 3). The CV of Chl. *a* concentrations in HSZ was similar to that in DLT (Table 3). The 12-month periodicity in Chl. *a* concentrations was significant for all reaches except for HSZ (Table 4, Figure 4). Moreover, HSZ had the highest mean Chl. *a* concentrations among the four reaches during the high flow period (Appendix A Figure A1); however, differences were not significant in this test (Kruskal–Wallis chi-squared = 5.092, df = 3, *p*-value = 0.165). During the low flow period, the mean Chl. *a* concentration in HSZ (7.17 mg m<sup>−</sup>2) was slightly higher than that in the tributary reach of DLT (5.59 mg m<sup>−</sup>2), and was lower than the other two mainstream reaches. These differences were also not significant (Appendix A Figure A1).

By using CCM, we found that water temperature had a significant influence on Chl. *a* dynamics in SNY, DLT, and JZG, but there was variation in the coupling of these variables through time and also the strength of this coupling (Figure 5; Appendix A Table A1; Appendix A Figure A2). The influence that water temperature had on Chl. *a* was tightly coupled in time except for in DLT where there was a 2-month delay in this coupling (Appendix A Table A1). By comparison, water temperature had the weakest and nonsignificant influence on Chl. *a* in HSZ as it had the lowest ρ value within the range of the 5th and 95th percentiles of the null model (Figure 5; Appendix A Table A1). In contrast to water temperature, discharge had a weak influence on Chl. *a* dynamics in all reaches, represented by low ρ values (Figure 5; Appendix A Table A1; Appendix A Figure A2). Furthermore, discharge did not affect Chl. *a* dynamics more significantly than that in the null model for all reaches (Figure 5).

**Figure 4.** Continuous wavelet power spectra showing the periodicity of benthic algal chlorophyll *a* concentrations (Chl. *a*) from July 2011 to June 2017 for each study reach. Further information on wavelet spectra are described in Figure 2. The original time series are demonstrated above each wavelet plot. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

**Figure 5.** *Cont*.

**Figure 5.** Convergent cross-mapping (CCM) between benthic chlorophyll *a* concentration (Chl. *a*) versus stream water temperature (WT) (left) and discharge (right) using time series from July 2011 to June 2017 for each study reach. The y-axis represents the Cross Map Skill (rho or Pearson's correlation coefficient) and the x-axis represents the library size (different lengths of data subsampled randomly from the time series) used for modeling. The time lags having the highest cross map skill in CCM were used in the final CCM for significance testing. Shaded areas are the 95% confidence intervals (5th to 95th percentiles) of 100 surrogate time series from the null model (light red area = 95% confidence intervals of the null model for the red line; light blue area = 95% confidence intervals of the null model for the blue line; overlap=light purple). rho improving with increasing library size indicates causality, and rho beyond the 95% confidence intervals of the null model implies significant causality. For JZG, the red line is beyond the light shaded area, indicating that WT has significant effects on Chl. *a*. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

#### **4. Discussion**

Differences in stream regime descriptors, including parameter magnitude, variability, frequency, timing, and predictability, among adjacent reaches can provide ecologically valuable information on how groundwater influences surface waters. We found that groundwater discharge can have a large influence on the temporal patterns of many stream regime metrics that describe variation in water temperature and discharge. Further, we found that changes in stream physical conditions due to groundwater discharge can alter the influence of these physical conditions for driving temporal variation in benthic algal biomass.

The most groundwater-influenced reach, HSZ, had a distinct air–water temperature relationship and a more stable thermal regime when compared with the other three reaches, confirming predictions and observations that this reach was substantially influenced by groundwater. In particular, HSZ had a lower maximum monthly mean water temperature than the adjacent higher elevation reach and a higher minimum monthly mean water temperature than the adjacent lower elevation reach. The mean water temperature in HSZ was also lower than that in higher elevation reach during the high flow period, and was higher than that in the lower elevation reach during the low flow period. Further, HSZ had the lowest frequency of months with temperatures < 3 ◦C among all reaches. These results indicate that stream water in HSZ was cooler in summer and warmer in winter compared to the other reaches. The smallest CV of water temperature and the shallowest slope for air–water temperature relationship in HSZ imply a slower heating and cooling rate (i.e., a temperature buffering effect) than the other reaches during the annual temperature cycle [20]. Such buffering effect has been commonly observed on daily temperature magnitude [54–56]. Our findings suggest that this phenomenon could also occur at longer time scales. The same timings for the maximum and minimum monthly mean water temperature across all reaches suggest that groundwater did not alter the timing of the metrics for water temperature in HSZ.

Although discharge increased with stream size, HSZ had the highest CV and the highest frequency for low-flow events. This result is related to the fact that groundwater was the sole source of flow for this reach when upstream reaches, which contain visible swallow holes, often ceased to flow between December to February. In fact, we found that mean discharge in HSZ (0.196 m<sup>3</sup> s−1) was a little smaller than the sum of flows in the two upstream reaches of SNY and DLT (0.225 m<sup>3</sup> s−1) in the whole study duration, and this difference was even more evident during the low flow period. It is likely that a proportion of downwelling flow had longer flow paths than the distance of the study reaches, longer subsurface residence times than we could detect in this study, or had lateral path(s) beyond the watershed boundary [57]. Therefore, although discharge magnitude increased longitudinally, karst topographic features likely induced a higher variability of flow dynamics in HSZ than in other reaches. Overall, the different regime metrics of discharge in HSZ compared with other reaches provides evidence that this reach is significantly influenced by groundwater discharge.

HSZ had a lower Chl. *a* magnitude and a higher low-concentration frequency than that in the upstream reach of SNY. This result seems unexpected because benthic algal biomass often increases with stream size and temperature [32]. Further, spring stream flow is generally rich in dissolved and organic nutrients that can promote downstream algal biomass [18,58]. One possible reason is that higher predation pressure by lotic grazers in HSZ reduced algal biomass because more stable temperature and flow conditions were also beneficial to grazer populations [59]. A comparison of the differences in nutrient concentrations between HSZ and nearby reaches and an evaluation of intensity of predation pressure and its influences on algal biomass in HSZ would be valuable topics for future research. Benthic algal biomass in HSZ did not display a significant annual cycle as that in other reaches. This phenomenon has previously been attributed to reduced temporal variation in physical conditions such as temperature, conductivity, and flow [60]. Our findings provide further support for this notion because water temperature and flow in HSZ had milder temporal regimes and their influences on temporal dynamics of Chl. *a* were weaker than reaches less influenced by groundwater. We found that stream discharge did not have a significant influence on the temporal dynamics of Chl. *a* at any study reach. A possible explanation for this is that stream discharge mainly affected benthic algal assemblages via indirect paths. Discharge was relatively low in most times during the study period and flow velocity may not have been high enough to reduce algal colonization or growth. Therefore, flow may have been most important for mediating the transfer of nutrients for autotrophic assimilation [31]. While floods might have negatively impacted lotic assemblages, the effect of these events may have been temporary and/or perhaps not detectible or important over log time periods, as in our study.

We used long-term data (six years) to gain a broad understanding of the physical and biotic influences of groundwater in our study stream. However, we were unable to adequately explore how groundwater–surface water exchange during contrasting hydrological conditions may influence important physical and biotic components of our stream system, due to the monthly sampling frequency of our study design. Further research with a more frequent and intensive survey effort would provide more valuable information on how the important physical and biotic components of streams that we assessed change in response to intra- and inter-annual variation in groundwater–surface water exchange.

In conclusion, water temperature in the most groundwater-influenced stream reach had the smallest variability (i.e., CV) and the lowest frequency for low-temperature months, when compared with other reaches. Groundwater inputs thus induced milder physical habitat conditions. These milder physical habitat conditions were associated with weaker temporal patterns for benthic algal biomass. That is, benthic algae biomass in this reach did not display a significant annual cycle as that in other reaches. Our findings demonstrate that groundwater inputs may alter regime metrics other than magnitude for physical components in downstream reach, and further induce substantial changes in temporal dynamics of important biotic components. Therefore, longer and more intensive studies will greatly improve our understanding of how groundwater effects stream ecosystems.

**Author Contributions:** Conceptualization, T.T. and Q.C.; Formal analysis, T.T.; Funding acquisition, T.T. and Q.C.; Investigation, S.G. and L.T.; Methodology, T.T.; Writing—original draft, T.T.; Writing—review & editing, T.L. and R.M.B.; and all authors reviewed the manuscript before submission.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 31470510), the National Key Research and Development Program of China (No. 2016YFC0503601, No. 2016YFC0502106), and the Special Project of Science and Technology Basic Work of Ministry of Science and Technology of China (No. 2014FY120200).

**Acknowledgments:** We thank Xiaoyu Dong, Tingting Sun, Ting Tang, Ze Ren, Fengzhi He, and Yang Li for assistance with fieldwork. We are grateful to Andrew Boulton and Naicheng Wu for valuable comments on the draft.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

**Table A1.** Convergent cross-mapping detecting the best cross-map skill (ρ) and time lag indicating the causal effects of stream water temperature (WT) or discharge on dynamics of benthic chlorophyll a concentration (Chl. *a*) from July 2011 to June 2017 for each reach. -: No causality is detected. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.


**Figure A1.** Plots demonstrating among-reach differences in the mean (± standard deviation) for water temperature, discharge, and chlorophyll a concentration (Chl. a) during the high-flow period (a, May to July) and the low-flow period (b, December to February), respectively. The χ<sup>2</sup> value and the *p*-value are the results of Kruskal–Wallis tests among the four reaches. The lower-case letters above whiskers are the results of pairwise comparisons, with the same letters indicating no difference at *p* ≤ 0.05. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

**Figure A2.** Extended convergent cross-mapping detecting causality and time lags between chlorophyll a concentrations (Chl. *a*) and water temperature (left) or discharge (right) using time series from July 2011 to June 2017 for each study reach. SNY = Sheng Nong Yaun; DLT = Da Long Tan; HSZ = Hong Shi Zi; JZG = Ji Zi Gou.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
