The Effects of Climate Variation and Anthropogenic Activity on Karst Spring Discharge Based on the Wavelet Coherence Analysis and the Multivariate Statistical

: This study focused on the impact of anthropogenic activity on magnitude, frequency, and minima of spring discharge. Niangziguan Springs (NS), China, was selected as an example, as its discharge is decreasing due to the combined effects of climate variation and human activity. For exploring the impact of human activity on the spring discharge from climate change, the spring discharges from 1959 to 2015 were divided into two periods: pre-development period (i.e., 1959–1980) and post-development period (i.e., 1981–2015). A polynomial regression model of the spring discharge was developed for the pre-development period. We deduced the model in the post-development period, compared the results with the observed spring discharge, and concluded that the climate variation and human activity caused 6.93% and 32.38% spring discharge decline, respectively. The relationships of spring discharge with Indian Summer Monsoon (ISM), East Asian Summer Monsoon (EASM), E1 Niño Southern Oscillation (ENSO), and Paciﬁc Decadal Oscillation (PDO) were analyzed by wavelet analysis during the two periods. The results illustrated that the monsoons (i.e., ISM and EASM) were dominated by climate factors that affect the NS discharge versus climate teleconnections (i.e., ENSO and PDO). According to different time scales, human activities have had an impact on the periodicity of NS discharge, which altered the periodicities of the spring discharge at inter-annual time scales, but the periodicities at intra-annual and annual time scales have remained the same between the two periods. Under the effects of human activity, the local parameter of non-stationary general extreme value (NSGEV) distribution varied with time. The predicted spring discharge minimum value is supposed to be 4.53 m 3 /s with a 95% conﬁdential interval with an upper boundary of 6.06 m 3 /s and a lower boundary of 2.80 m 3 /s in 2020. The results of this study would beneﬁt the management of spring discharge and water resources.


Introduction
Carbonate rocks are easily dissolved by groundwater, and an abundant amount of voids, fractures, and conduits are formed. Such an aquifer is called a karst aquifer. Around one-tenth of the Earth's surface is covered by carbonate rocks, which serve as a crucial source of water, providing hydration to almost 25% of the global population [1]. Karst aquifers are, however, particularly vulnerable to environmental changes [2]. During the last few decades, climate variation and intensive groundwater development have led to many environmental problems for karst aquifers, such as groundwater level decrease, pumping costs increase, natural spring dry-up, land subsidence, and ecosystem damage [3]. These problems have promoted karst groundwater research in many places throughout (PDO) in two different time periods. By comparing the results between the two periods, we assess the impacts of climate variation and human activity on spring discharge in terms of magnitude, frequency, and minima. This study demonstrates the impacts of climate variation and human activity on spring discharge at multiple scales, and its findings have significant implications for the sustainable utilization and environmental protection of karst water resources.

Site Description
The NS basin is situated along the Mianhe riverbank, spanning 7 km in the Mianhe Valley of Eastern Shanxi Province (Figure 1). The NS complex discharge from Ordovician karstic aquifers, which are overlaid by Carboniferous sand shales, Permian and Triassic detrital formations, and Quaternary deposits, upward [25]. Karst groundwater converges from the north, west, and south toward Mianhe Valley in the east. In this location, groundwater accumulates in the low-permeable strata of Cambrian dolomite, which intersects with the surface to create NS [26,27]. Because the groundwater in the NS basin is separated with those in adjacent regions by an impermeable boundary, a low-permeable fault boundary, tectonic divide, groundwater divide, and surface water divide, the groundwater in the basin has no water exchange with neighboring regions [26]. Precipitation is the principal means of recharging karst aquifers, and spring discharge and groundwater pumping are the two major ways of groundwater consumption in the NS basin. (EASM), Indian Summer Monsoon (ISM), El Niño-Southern Oscillation (ENSO), and Pacific Decadal Oscillation (PDO) in two different time periods. By comparing the results between the two periods, we assess the impacts of climate variation and human activity on spring discharge in terms of magnitude, frequency, and minima. This study demonstrates the impacts of climate variation and human activity on spring discharge at multiple scales, and its findings have significant implications for the sustainable utilization and environmental protection of karst water resources.

Site Description
The NS basin is situated along the Mianhe riverbank, spanning 7 km in the Mianhe Valley of Eastern Shanxi Province (Figure 1). The NS complex discharge from Ordovician karstic aquifers, which are overlaid by Carboniferous sand shales, Permian and Triassic detrital formations, and Quaternary deposits, upward [25]. Karst groundwater converges from the north, west, and south toward Mianhe Valley in the east. In this location, groundwater accumulates in the low-permeable strata of Cambrian dolomite, which intersects with the surface to create NS [26,27]. Because the groundwater in the NS basin is separated with those in adjacent regions by an impermeable boundary, a low-permeable fault boundary, tectonic divide, groundwater divide, and surface water divide, the groundwater in the basin has no water exchange with neighboring regions [26]. Precipitation is the principal means of recharging karst aquifers, and spring discharge and groundwater pumping are the two major ways of groundwater consumption in the NS basin.  According to records from 1959 to 2015, NS has been identified as the most extensive karst springs in northern China, boasting an average annual discharge of 9.81 m 3 /s ( Figure 2). The catchment area of the springs encompasses 7394 km 2 , which includes Yangquan City, as well as the counties of Pingding, Heshun, Zuoquan, Xiyang, Yuxian, Yangquan, and Shouyang, and serves as the main water source for the NS (Figure 1). The NS basin is situated in a semiarid continental monsoon climate within the warm temperate zone. Thus, the area is hot and rainy in summer, and cold and dry in winter. There are significant differences in precipitation the whole year, with the main rainfall concentrated from July to September in summer (the maximum annual precipitation is 843.85 mm). On the contrary, there is very little rainfall in winter, with a minimum annual precipitation of 292.57 mm. Based on records from 1959 to 2015, the average annual precipitation is 490.5 mm (Figure 2). According to records from 1959 to 2015, NS has been identified as the most extensive karst springs in northern China, boasting an average annual discharge of 9.81 m 3 /s ( Figure  2). The catchment area of the springs encompasses 7394 km 2 , which includes Yangquan City, as well as the counties of Pingding, Heshun, Zuoquan, Xiyang, Yuxian, Yangquan, and Shouyang, and serves as the main water source for the NS (Figure 1). The NS basin is situated in a semiarid continental monsoon climate within the warm temperate zone. Thus, the area is hot and rainy in summer, and cold and dry in winter. There are significant differences in precipitation the whole year, with the main rainfall concentrated from July to September in summer (the maximum annual precipitation is 843.85 mm). On the contrary, there is very little rainfall in winter, with a minimum annual precipitation of 292.57 mm. Based on records from 1959 to 2015, the average annual precipitation is 490.5 mm ( Figure 2). Since the 1980s, the rapid development and exploitation of karst water resources in the NS basin have been identified as the cause of the decrease in NS discharge [28,29]. The observed increase in karst water utilization can be primarily attributed to four major human activities. The foremost factor behind the significant rise in karst water usage was dewatering activities associated with coal mining operations in the region. The impact of these activities on the spring discharge has been quantified at 1.02 m 3 /s [30]. The second significant human activity was the construction of anti-seepage channels and dams. The creation of anti-seepage channels served to prevent surface water from infiltrating into the subsurface, while the construction of dams helped to store extensive amounts of water Since the 1980s, the rapid development and exploitation of karst water resources in the NS basin have been identified as the cause of the decrease in NS discharge [28,29]. The observed increase in karst water utilization can be primarily attributed to four major human activities. The foremost factor behind the significant rise in karst water usage was dewatering activities associated with coal mining operations in the region. The impact of these activities on the spring discharge has been quantified at 1.02 m 3 /s [30]. The second significant human activity was the construction of anti-seepage channels and dams. The creation of anti-seepage channels served to prevent surface water from infiltrating into the subsurface, while the construction of dams helped to store extensive amounts of water from the catchment; as a result, this has led to a reduction in river flow and a decrease in groundwater recharge within the exposed karst valley. The impact of these activities on karst water recharges has been estimated at a decline of 0.84 m 3 /s [31]. The third significant human activity attributed to the increase in karst water utilization was the extraction and exploitation of karst groundwater resources. The primary impact of this activity has been in decline in spring discharge estimated at 0.46 m 3 /s [28]. The fourth primary human activity leading to an increase in karst water usage was the exploitation of unconsolidated sediment aquifers situated in the upper regions. The impact of this activity has resulted in a decline in spring discharge estimated at 0.06 m 3 /s [28]. Previous studies have presented evidence supporting the segmentation of spring discharge data.

Data Sources
From January 1959 to December 2015, Niangziguan Spring (NS) discharge monthly data were collected from the Niangziguan gauge station. During the same period, monthly precipitation data were obtained from six meteorological stations located in the NS basin. Additionally, monthly temperature data were collected from meteorological stations situated within the Yangquan area for the same duration ( Figure 2).
The monthly average data for the EASM were gathered from the National Oceanic and Atmospheric Administration's (NOAA) website (http://www.cpc.ncep.noaa.gov (accessed on 1 December 2022)). In addition, the Monsoon Monitoring website (http://apdrc.soest. hawaii.edu (accessed on 1 December 2022)) provides monthly average data for monsoon indices of ISM. Moreover, the monthly data for the teleconnection patterns of ENSO were obtained from the Climate Prediction Center of the National Weather Service website, whereas data for the Pacific Decadal Oscillation (PDO) are available from the Joint Institute of the Study of the Atmosphere and Ocean, on the University of Washington (JISAO) website (http://jisao.washington.edu (accessed on 1 August 2022)).

Methods
In order to analyze these datasets, we used principal component analysis, cluster analysis and the significance test, continuous wavelet transform (CWT), wavelet coherence, global coherence, and non-stationary general extreme value (NSGEV) distribution. Brief descriptions of these methods and the reasons behind employing them are discussed below.

Principal Component Analysis (PCA)
Seven precipitation stations were distributed within the NS basin. Thus, the precipitation data from seven stations were subjected to principal component analysis, and the first principal component derived from original data was employed to represent precipitation in the NS Basin. It was the precipitation dataset for the rest of analysis. A brief description of the principle of the PCA is given below.
PCA is a statistical method that requires an orthogonal transformation to convert a dataset of inter-connected variables into a set of linearly uncorrelated variable values [32]. In general, the statistical procedure of PCA can reveal the internal structure of the original datasets. PCA frequently involves converting a high-dimensional dataset into a lowdimensional dataset. This is achieved by selecting only the first few principal components, resulting in a reduction in the dimensionality of the transformed data [33].

Cluster Analysis and the Significance Test
The NS basin was a relatively underdeveloped region in China until the 1980s. In the last 30 years, it has gone through a notable transformation and developed into one of the leading industrialized centers in the nation, particularly in the areas of coal mining, power generation, and metallurgy. Thus, the monthly spring discharge data from January 1959 to December 2015 were analyzed by cluster analysis to examine the impact of human disturbance on spring discharge, and the rationality of the section was analyzed by a Wilcoxon rank-sum test. Briefly, cluster analysis aims to categorize a collection of observational samples based on their similarities, in such a way that samples within the same group are more comparable to each other than to those outside the group [34]. In other words, the purpose of cluster analysis is to divide samples into groups, where samples of the same group have some properties in common. Initially, each sample is regarded as a group, and then the distances between groups are calculated by a median method. The pairs of the samples which have the minimum distance are merged into a new group, and then the distance between the new group and other groups are recalculated. This process continues recursively until no samples change groups.
The Wilcoxon rank-sum test is a nonparametric statistical test utilized to examine the null hypothesis, and it is capable of determining the significance of the discrepancy between two sample series. After combining the two sample series, the observation series are ranked in ascending order and assigned their respective ranks. In the Wilcoxon rank-sum test, the rank-sums of observations in two sample sets are computed, and significant levels between the two sets are assessed with a hypothesis test [35]. The Wilcoxon rank-sum test relies exclusively on the ordering of the observed values between the two sample sets [36,37]. Unlike the t-test, it does not necessitate the data to follow a normal distribution. However, it is almost as effective as the t-test when the data distribution is normal.

Continuous Wavelet Transform (CWT)
The continuous wavelet transform (CWT) is a highly favored tool amongst timefrequency transformations [38]. It is a mathematical operation that convolves the input data sequence with the set of functions created from the "mother wavelet". This analysis method is commonly used in wavelet analysis, which decomposes signals into a series of wavelet functions, characterized by irregular and asymmetrical time scales. A wavelet is a localized function, having a distinct number of oscillations that persist over a specific time span before rapidly fading towards zero.
The unique characteristics of the CWT sets it apart from the trigonometric functions employed in conventional Fourier analysis [39]. Wavelet transforms can be utilized to analyze the temporal and frequency variations in both the NS discharge and climate indices [25,40]. CWT method often exhibits edge artifacts since the wavelet used is not entirely localized in time. To address this issue, the Cone of Influence (COI) is usually introduced. COI is a critical region in the wavelet spectrum where the edge effects of the transform cannot be disregarded. By incorporating COI, the effects of edge artifacts can be mitigated, resulting in more accurate and reliable wavelet analysis. Wavelet analysis was used to detect the significant periods and regions of climate indices in the time-frequency domain, and provided the coefficients for calculating the wavelet coherence coefficients.

Wavelet Coherence
Wavelet coherence is utilized to inspect the correlation between climate indices and spring discharge in the time-frequency domain, aiming to identify the periods and time scales of high correlation between the two series [40]. It is defined the coefficient of wavelet coherence as: where a is a scale parameter in the frequency domain; τ denotes a positional parameter within the time domain. W x (a, τ) is the wavelet power that reflects essential characteristics of a time series x t , while W y (a, τ) reflects essential characteristics of a time series y t . W xy (a, τ) denotes the cross-wavelet power, which can be represented as the crosscorrelation coefficients. S represents a smoothing operator. R 2 (a, τ) as wavelet coherence measurement typically ranges between 0 to 1, indicating the strength of correlation in the time-frequency domain of the analyzed signal. This value represents the localized correlation coefficient obtained through wavelet analysis. To derive this coefficient, a smoothing operator referred to as S is defined: where S scale denotes the smoothing conducted along the wavelet scale axis and S time represents smoothing in the time scale. During practical implementation, both S scale and S time convolutions are computed discretely. Consequently, the computation necessitates numerical determination of normalization coefficients. Monte Carlo simulation method is utilized to determine the statistical significance of the wavelet coherence. The computation of wavelet coherence is carried out for every surrogate pair, and the estimation of the significance level for each scale takes into account the values that are situated beyond the COI. Empirical findings suggest that using 12 scales per octave is adequate. In accordance with Torrence and Compo's methodology, our analysis considers a 5% significance level against red noise [40].

Global Coherence
The determination of the global wavelet coherence coefficient for a specific scale a can be achieved by taking the average of the wavelet coherence coefficients across time: where n represents the quantity of points within the time scale [40,41]. The global coherence coefficient is a statistical tool that can measure the correlation between two time series, with the added benefit of being able to measure on multiple scales. This is an effective way for analyzing data patterns and periodicities without being affected by time [40,42].

Non-Stationary General Extreme Value (NSGEV) Distribution
To assess the risk of the NS drying up, we analyzed the minimum spring discharge data using the NSGEV distribution [43]. NSGEV distribution is used to fit the non-stationary time series constituted by maximum values whose structure is described as follows: and its probability distribution function is: where µ t , σ t , and ξ t represent the location, scale, and shape parameters, respectively, and In general, µ t and σ t are polynomial functions [43,44]. For a minimum sequence, X t , its opposite number sequence − X t becomes a maximum sequence, which also satisfies expression (4). A common approach for modeling it is by with distribution function where µ t = −µ t , σ t = σ t , and ξ t = ξ t , and {x :

PCA of Precipitation Dataset
As mentioned before, the precipitation data for the study area were gathered from six meteorological stations situated within the NS basin ( Figure 1). PCA was used to examine the correlations between six precipitation sequences and extract the maximum value of the precipitation characteristics. Table 1 presents the results of PCA conducted on a dataset consisting of six precipitation sequences. The analysis yielded a first principal component (Y1) with a variance proportion of 90.29%, exceeding the minimum threshold value of 80%. For this reason, Y1 was used as a representation of precipitation conditions in the NS basin.

The Spring Discharge Segment and the Significant Test
Using cluster analysis, we classified the NS discharge data, and obtained a hierarchical clustering which was a set of nested clusters that were organized as a binary tree. Each node of the tree represented a separation year of its sub-clusters. The node of the first level was 1980, which indicated that the NS discharge could be divided into two groups at 1980. The region's economic growth has led to significant water resource consumption. Prior to 1980, surface water sources served as the primary water supply for local stations in the NS basin. However, the rapid economic development since the 1980s has primarily relied on groundwater resources. Therefore, the cluster of the NS discharge at 1980 is reasonable, and the data of the NS discharge can be divided into two groups and the timeline for groundwater in the NS basin can be divided into two distinct periods. The pre-development period spanned from 1959 to 1980, where groundwater was mainly influenced by climatic variations, and the impact of human activities was considered negligible. The post-development period lasted from 1981 to 2015, during which both climatic and anthropogenic factors were identified as significant contributors to changes in groundwater quality. We utilized the Wilcoxon rank-sum test to assess the importance of the spring discharge segment in 1980 [36]. The p-values obtained for the yearly average spring discharge were 2.2 × 10 −16 . This indicates that the distributions of spring discharge were significantly distinct, comparing the periods before and after the development. The value was far below the threshold of 0.05.

Evaluation of the Reduction in Spring Discharge Attributable to Human Activity
In the NS basin, the annual average precipitation was 569 mm in the pre-development period, and 500 mm in the post-development, showing a 12.03% of decrease for precipitation. Similarly, during the pre-development period, the average yearly temperature was 10.9 • C. In contrast, it was 11.4 • C in the post-development time, constituting an increase of 4.59%.
The annual average spring discharge during the pre-development period was measured at 12.54 m 3 /s, which significantly reduced to 7.61 m 3 /s during the post-development period, culminating in a decline of 39.31%. In order to separate the contribution to the decline of the spring discharge due to human activity from that due to climate change, we analyzed the relations of the spring discharge with precipitation and temperature of the pre-development period, using a polynomial regression [45]. The analysis showed that temperature was not significantly related with the spring discharge. However, the relations between spring discharge and precipitation in pre-development were: Q(t) = 8.32 + 4.68 × 10 −3 P t + 7.01 × 10 −9 P 3 t−1 (8) where Q(t) denotes the annual average spring discharge at the time of t (m 3 /s); P t is annual precipitation at the time of t (mm); P t−1 means annual precipitation at the time of t − 1 (mm). We then estimated the spring discharge under effects of climate change in postdevelopment period by using Equation (8). The annual average spring discharge only under effects of climate change in post-development period was 11.67 m 3 /s. Thereby, the contribution of climate change to the spring discharge decline was likely to be 0.87 m 3 /s (i.e., 12.54 m 3 /s − 11.67 m 3 /s = 0.87 m 3 /s), and that of human activity to be 4.06 m 3 /s (i.e., 11.67 m 3 /s − 7.61 m 3 /s = 4.06 m 3 /s). Such a result indicates that the climate variation and human activity contribute 6.93% and 32.38% to spring discharge decline, respectively.

CWT of the Climate Phenomenon, Precipitation and the Spring Discharge
CWT technique was implemented to detect the primary modes of oscillation within time series of monthly climate influences, such as ISM, EASM, ENSO, PDO, precipitation, and spring discharge. It is worth mentioning that a significant level of 0.05 was used for periodicities greater than or equal to one year, while a more stringent significant level of 0.01 was applied to periodicities less than one year. The selection of the significant level was aimed at excluding white noise of high frequencies. Figure 3a showed that ISM exhibited a continuous periodicity of 1 year during the pre-and post-development periods. Figure 3b illustrated that EASM mainly possessed continuous periodicity of 1 year during the both periods. Figure 3c showed that the periodicities of ENSO mainly concentrated at the time scales of 1 to 5 years in the preand post-development periods. Figure 3d illustrated the periodicities of PDO mainly focused on 0.5 years in pre-development period, and 0.5 to 1 year in the post-development period. Thus, for the monsoons, both ISM and EASM exhibited an obvious periodicity of 1 year, and for the climate teleconnections, ENSO has periodicities of 1 to 5 years, and PDO possesses periodicities of 0.5 to 1 year. Moreover, the ISM, EASM, and ENSO keep the same periodicities between the pre-and post-development periods, while PDO is a little different ( Table 2). Because most of the climate indices keep the same periodicities, one can detect the spring discharge changes due to human activity by correlations of the climate indices with spring discharge between the two periods. Table 2. The periodicity of indices in the pre-development and post-development period by wavelet analysis (year).

Indices
Pre-Development Post-Development  (Table 2). In karst hydrological processes, precipitation infiltrates to the subsurface, becoming groundwater, which flows through voids, fractures, and conduits, and then reaches outlets discharging as springs. In the process, the precipitation signals are attenuated due to resistance and storage of the karst aquifer. Only the strong precipitation signals can propagate through the karst aquifer and emerge as spring discharge. For this reason, the CWT coefficients of the spring discharge were weaker than those for precipitation. Even though the spring discharge signals were weak, the correlations of climate indices and precipitation with spring discharge were still significant. In other words, one can identify the spring discharge changes due to human activity by comparing the wavelet coherence coefficients of spring discharge with climate indices between pre-and post-development periods.   storage of the karst aquifer. Only the strong precipitation signals can propagate through the karst aquifer and emerge as spring discharge. For this reason, the CWT coefficients of the spring discharge were weaker than those for precipitation. Even though the spring discharge signals were weak, the correlations of climate indices and precipitation with spring discharge were still significant. In other words, one can identify the spring discharge changes due to human activity by comparing the wavelet coherence coefficients of spring discharge with climate indices between pre-and post-development periods.  Figure 5a shows the continuous periodicity of 1 year and the intermittent periodicity of 0.5 and 2 years exhibited between precipitation and ISM in the pre-development period. The continuous periodicity of 1 year and intermittent periodicity of 0.5 and 2.5 years is exhibited between precipitation and ISM in the post-development period (Figure 5b). Moreover, the global coherence coefficients of precipitation with ISM at periodicities of 0.5 and 1 year were very close between the two periods, respectively (Figure 5e). However, the differences of the global coherence coefficients of precipitation with ISM were observed at inter-annual time scales between the two periods ( Figure 5e). The results  Figure 5a shows the continuous periodicity of 1 year and the intermittent periodicity of 0.5 and 2 years exhibited between precipitation and ISM in the pre-development period. The continuous periodicity of 1 year and intermittent periodicity of 0.5 and 2.5 years is exhibited between precipitation and ISM in the post-development period (Figure 5b). Moreover, the global coherence coefficients of precipitation with ISM at periodicities of 0.5 and 1 year were very close between the two periods, respectively (Figure 5e). However, the differences of the global coherence coefficients of precipitation with ISM were observed at inter-annual time scales between the two periods ( Figure 5e). The results indicated that ISM affects precipitation stably at intra-annual and annual time scales, but impacts precipitation at inter-annual time scale unstably.

Wavelet Coherence of Precipitation and Spring Discharge with ISM
Similarly, the correlation between spring discharge and ISM is exhibited during the pre-development period in Figure 5c, particularly in terms of intermittent periodicity around 0.5 and 1 year. On the other hand, the intermittent periodicity around 0.5, 1, 2.5, and 5 years is exhibited in the post-development period (Figure 5d). Further, the global coherence coefficients of the spring discharge with ISM around periodicities of 0.5 and 1 year were observed in both pre-and post-development periods, which means that human activity could not alter the resonant periodicities of the spring discharge with ISM at intraannual and annual time scales (Figure 5f). However, the inconsistencies of global coherence coefficients of the spring discharge with ISM at inter-annual time scales were observed between the two periods. This finding indicates that human activity caused a variation in the resonance patterns of the spring discharge with the ISM. In addition, the global coherence coefficients of precipitation with ISM were much larger than that of the spring discharge with ISM (Figure 5f). The results suggest that when precipitation infiltrated into the subsurface and propagated through the aquifer, its signals were attenuated. annual and annual time scales (Figure 5f). However, the inconsistencies of global coherence coefficients of the spring discharge with ISM at inter-annual time scales were observed between the two periods. This finding indicates that human activity caused a variation in the resonance patterns of the spring discharge with the ISM. In addition, the global coherence coefficients of precipitation with ISM were much larger than that of the spring discharge with ISM (Figure 5f). The results suggest that when precipitation infiltrated into the subsurface and propagated through the aquifer, its signals were attenuated.   Figure 6a displays the continuous periodicity of 1 year and the intermittent periodicity of 0.5 and 8 years exhibited between precipitation and EASM in the pre-development period. The continuous periodicity of 1 year and the intermittent periodicity of 0.5, 2, 3, and 8 years is exhibited between precipitation and EASM in the period after development (Figure 6b). Moreover, the coherence coefficients at 0.5-and 1-year periodicities between EASM and global precipitation were highly similar in both periods (Figure 6e). Meanwhile, during the two observed periods, variations in global coherence coefficients of precipitation and EASM were noted at inter-annual time scales. The observations implied that the impact of EASM on precipitation was consistent on an annual and intra-annual basis; its impact becomes unpredictable at an inter-annual time scale. and 8 years is exhibited between precipitation and EASM in the period after development (Figure 6b). Moreover, the coherence coefficients at 0.5-and 1-year periodicities between EASM and global precipitation were highly similar in both periods (Figure 6e). Meanwhile, during the two observed periods, variations in global coherence coefficients of precipitation and EASM were noted at inter-annual time scales. The observations implied that the impact of EASM on precipitation was consistent on an annual and intra-annual basis; its impact becomes unpredictable at an inter-annual time scale. Similarly, Figure 6c illustrates the intermittent periodicity around 0.4 and 1 year exhibited between spring discharge and EASM during the pre-development period. On the other hand, the intermittent periodicity around 0.4, 1, 2.4, and 7 years is exhibited between spring discharge and EASM in the post-development period (Figure 6d). The coherence coefficients for the spring discharge with EASM showed consistent global values in the Similarly, Figure 6c illustrates the intermittent periodicity around 0.4 and 1 year exhibited between spring discharge and EASM during the pre-development period. On the other hand, the intermittent periodicity around 0.4, 1, 2.4, and 7 years is exhibited between spring discharge and EASM in the post-development period (Figure 6d). The coherence coefficients for the spring discharge with EASM showed consistent global values in the two periods, specifically around the 0.4-year and 1-year periodicities, which suggested that human activity could not alter the resonant periodicities of the spring discharge with EASM at intra-annual and annual time scales (Figure 6f). However, the inconsistencies of the global coherence coefficients of the spring discharge with EASM were observed at interannual periodicities between two periods, which revealed that the resonant periodicities of the spring discharge with EASM were altered by human activity. ter-annual periodicities between two periods, which revealed that the resonant periodici-ties of the spring discharge with EASM were altered by human activity. Figure 7a shows the intermittent periodicities around 0.4, 0.8, and 2 years exhibited between precipitation and ENSO in the pre-development period. The intermittent periodicity around 0.4, 0.8, and 3 years is exhibited between precipitation and ENSO in the post-development period (Figure 7b). Moreover, the global coherence coefficients of precipitation with ENSO around periodicities of 0.4, 0.8, and 2-3 years were close between the two periods, respectively (Figure 7e). The results indicated that ENSO mainly affected precipitation in the NS basin at intra-annual and inter-annual time scales.  Similarly, Figure 7c,d illustrates the intermittent periodicities around 0.4, 0.8, and 3 years exhibited between the spring discharge and ENSO in the two periods. Moreover, the global coherence coefficients around periodicities of 0.4, 0.8, and 3 years were observed in the two periods, respectively (Figure 7f). Through comparing the global coherence coefficients of the spring discharge with ENSO between pre-and post-development, one can conclude that the human activity did not alter the resonant periodicities of ENSO with the spring discharge. Figure 8a,b shows the intermittent periodicities around 0.5 and 1 year exhibited between precipitation and PDO in the two periods. Although the periodicity of 5 years was still observed between the two periods, it was out of the COI in pre-development, and did not pass the test of 0.05 significance level in post-development. Thus, the resonant periodicity of 5 years was negligible. Figure 8e illustrated that the global coherence coefficients of precipitation with PDO around periodicities of 0.5 and 1 year were consistent between the two periods. However, the coefficient values were small with large differences, which indicated that the effects of PDO on precipitation were weak in the NS basin. efficients of the spring discharge with ENSO between pre-and post-development, one can conclude that the human activity did not alter the resonant periodicities of ENSO with the spring discharge. Figure 8a,b shows the intermittent periodicities around 0.5 and 1 year exhibited between precipitation and PDO in the two periods. Although the periodicity of 5 years was still observed between the two periods, it was out of the COI in pre-development, and did not pass the test of 0.05 significance level in post-development. Thus, the resonant periodicity of 5 years was negligible. Figure 8e illustrated that the global coherence coefficients of precipitation with PDO around periodicities of 0.5 and 1 year were consistent between the two periods. However, the coefficient values were small with large differences, which indicated that the effects of PDO on precipitation were weak in the NS basin.  As shown in Figure 8c, in the pre-development period, there were intermittent periodicities at 0.5-, 1-, and 2.5-year intervals in the correlations between spring discharge and PDO. However, in the post-development period, there were intermittent periodicities at 1-and 2-year intervals in the correlations between spring discharge and PDO (Figure 8d). Figure 8f illustrates that the global coherence coefficients of the spring discharge with PDO around periodicities of 1, and 2 years were consistent between the two periods. However, the coefficient values were small with large differences, which demonstrated that the effects of PDO on the spring discharge were weak.

Wavelet Coherence of Precipitation and Spring Discharge with PDO
By comparing the global coherence coefficients of both monsoons (ISM and EASM) and climate teleconnections (ENSO and PDO) with precipitation, one can find that the former was much larger than latter. The results emphasized that the monsoons were dominate climate factors that affected precipitation in the NS basin. In a similar way, the monsoons also were dominate climate factors that affected the spring discharge.

NSGEV for the Spring Discharge Minima
To investigate the spring discharge minima distribution during the pre-development and post-development periods, NSGEV models were set up for such two periods, respectively. The parameters of NSGEV distribution were identified by the maximum likelihood. NSGEV distribution for the pre-development period was where Q 1 represents the spring discharge minima during the pre-development period; G(Q 1 ) represents the spring discharge minima distribution during the pre-development period. The NSGEV model for the post-development period was set up by using the spring discharge minima from 1981 to 2012. The remaining spring discharge minimum values from 2013 to 2015 were used to validated the model (Table 3).
where Q 2 represents the spring discharge minima during the pre-development period; G(Q 2 ) represents the spring discharge minima distribution during the post-development period. Equation (9) showed that local parameter µ t (i.e., 11.751), scale parameter σ t (i.e., 1.460), and shape parameter ξ t (i.e., −0.525) were constant during the pre-development period. In other words, the spring discharge minima were stationary under natural conditions. Equation (10) showed that local parameter µ t (i.e., 8.662 − 0.095t) varied with time, and that scale parameter σ t (i.e., 0.805) and shape parameter ξ t (i.e., −0.212) were constant during the post-development period. The downtrend of the local parameter indicated that human activity caused the spring discharge minima to decline.
Using Equation (10), we predicted the spring discharge minima from 2013 to 2020 (Table 3). Table 3 shows that the spring discharge minima (i.e., from 2013 to 2015) used for model validation would fall in the 95% confidential interval. Thus, the model was accurate, and the results were confidential. Results also showed that on average, the spring discharge minima was supposed to be 4.53 m 3 /s in 2020, and its upper boundary of 95% confidence was supposed to be 6.07 m 3 /s, and its lower boundary of 95% confidence was supposed to be 2.80 m 3 /s.

Conclusions
Since the 1980s, the rapid development and exploitation of karst water resources in the Niangziguan basin have led to a decline in spring discharge. In order to better analyze the impact of human activities and climate variation on spring discharge, this study distinguished the monthly spring discharge data from January 1959 to December 2015 through cluster analysis. From the investigation, we conclude that the climate variation accounted for 6.93% of the NS discharge decline, while human activity caused 32.38% of the NS discharge decline on average since 1981.
Compared with the results of the wavelet analysis during the two periods, it can be seen that the monsoons (ISM and EASM) were dominate climate factors that affected the NS discharge versus climate teleconnections (ENSO and PDO). Most of the climate indices (i.e., ISM, EASM, and ENSO) keep the same periodicities. A comparison of wavelet coherence coefficients of spring discharge with climate indices during pre-and post-development periods can help to identify the impact of human activity on spring discharge. Human activity affected the periodicity of NS discharge depending on the time scales. That is to say, human activity could alter the resonant periodicities of the spring discharge with climate indices at inter-annual time scales, but the resonant periodicities at intra-annual and annual time scale remained the same between the two periods.
Human activity caused the fact that the local parameter of (non-stationary general extreme value) NSGEV distribution varied with time in a downward trend. The spring discharge minima was supposed to be 4.53 m 3 /s with a 95% confidential interval with an upper boundary of 6.06 m 3 /s and a lower boundary of 2.80 m 3 /s in 2020.
A karst spring is the natural discharge point of an aquifer, and its discharge variation reflects the regional change of groundwater level. In general, the methods developed in this study can be applied to many other karst springs in northern China. This study has taken a solid step towards utilizing meteorological changes to predict groundwater resources variation. However, a karst aquifer is generally highly uneven. The impact of climate change and human activities on spring discharge often lags behind, and we will explore this topic in more depth in the future.