Delineation of Hydraulic Flow Regime Areas Based on the Statistical Analysis of Semicentennial Shallow Groundwater Table Time Series

: Shallow groundwater acts as an important source of water for the ecosystem, agriculture, drinking water supply, etc.; it is, however, among those water resources most sensitive to climate change, and especially to aridiﬁcation. In the present study, the delineation of regional recharge and discharge zones of the Danube–Tisza Interﬂuve (Hungary, 8000 km 2 ) is presented via the combination of multivariate time series and geomathematical methods to explore the subregions most sensitive to dewatering. The shallow groundwater level time series of 190 wells, covering a semicentennial period (1961 to 2010), were grouped into three validated clusters representing characteristically di ﬀ erent subregions. Then, the subregions’ means and individual shallow groundwater level time series were investigated for long-term trends and compared with local meteorological variability (precipitation, evapotranspiration, etc.) to determine their regime characteristics. As a result, shallow recharge and discharge zones, a gravity-driven ﬂow system, and the discharge zone of a deeper, overpressured ﬂow system could be discerned with distinctive long-term changes in water levels. The semicentennial trends in shallow groundwater levels were signiﬁcant (p < 0.05) in the recharge ( − 0.042 m y − 1 ) and in the overpressured discharge zone (0.009 m y − 1 ), and insigniﬁcant in the rest of the area ( − 0.005 m yr − 1 ). The present results concur with previous ﬁndings from the area but provide a statistically sound and reproducible delineation of the regime areas on a much ﬁner scale than before. With the determination of the di ﬀ erent climatic processes driving the semicentennial trends prevailing in the shallow groundwater, the high vulnerability of the recharge zone is underlined, while the outlined overpressured ﬂow system seems to act independently from semicentennial precipitation trends. This study provides a more in-depth picture of the long-term changes in shallow groundwater and its drivers in of one of the most important agricultural areas in Hungary. It outlines, in a generally applicable way, the most vulnerable subareas for irrigation relaying on shallow groundwater extraction. In addition, the results can help adaptation-strategy decision makers to initiate a more e ﬀ ective and area-focused intervention in the case of the predicted negative trends for vulnerable recharge areas under various climate change scenarios.


Introduction
Shallow groundwater is a dominant factor in sustaining river baseflow and global ecosystem functions [1]. Its quality and quantity are, however, being increasingly threatened by urbanization, Water 2020, 12 industrialization, and climate change [2,3]. In the past 50 years, groundwater extraction has reached a level of 35% of water withdrawn globally and has become an important source of fresh water, worldwide, for industrial (27%), domestic (36%), and agricultural (42%) purposes [4]. Cumulative global groundwater depletion is estimated to be between 4000 and 20,000 km 3 , corresponding to a sea level rise of 12 to 55 mm since 1950, depending on the estimation method. Furthermore, studies from around the world have shown that the increasing mean temperatures and changing rainfall patterns associated with climate change have a significant negative impact on groundwater recharge (i.e., long-term average renewable groundwater resources) and accessibility [2,[5][6][7]. As a result of these factors, in many countries, the depth of the groundwater table has started to decrease in the last few decades [8] (e.g., the United States, Australia [9], and China [10,11]). Groundwater recharge is also likely to be reduced in central and eastern Europe, with an even greater reduction in valleys and lowlands [3] (e.g., in the Hungarian steppes) [12]. Shallow groundwater (SGW) is an important source for irrigation and food production in many countries, and also serves as an important secondary source of water supply during natural droughts [13]. Furthermore, the state of the SGW table plays an important role in groundwater-dependent ecosystems [14], and thanks to capillary action, SGW, and therefore its level, also has an effect on surface vegetation and yield of agricultural crops [15,16]. Consequently, any decline in the water table, due to changes in precipitation patterns and rises in temperature, can reduce well yield and increase pumping costs. This could well have a serious economic impact on the livelihood of farmers in the regions where groundwater is used as the major source of irrigation [6,17].
Long term changes in SGW levels induced by climate change depend mainly on: • the hydraulic regime characteristics of the given area [18][19][20][21], as recharge areas are more sensitive to changes in local meteorological conditions than discharge ones; • the scale of the subsurface flow, while local flow systems are more sensitive to changes in infiltration conditions and changes in local climate, regional flow systems operating over a larger spatial scale (subregional and regional) are more resistant to such phenomena [22]; • the characteristic hydrogeological parameters of the aquifer, e.g., hydraulic conductivity, specific yield for unconfined aquifers, and specific storage for confined ones [14,[20][21][22][23].
In the short term (at the level of interannual variability), changes in SGW levels are driven primarily by evapotranspiration and the interannual distribution of precipitation [23]. Thus, finding significant correlations between precipitation and SGW level variations can help assess aquifer vulnerability [24]. The amount and the annual distribution of precipitation and evapotranspiration is mostly determined by the water balance of the recharge areas, while the SGW in the discharge areas can be replenished from subsurface recharge originating from the surface recharge zones, buffering the negative effect of precipitation loss in certain years [20,22].
Considering these factors together assists in the determination of the different flow regime areas (recharge and discharge), enabling the delineation of those areas most sensitive to climate change or severe anthropogenic impacts [14]. In the case of gravitational flow systems, the SGW table usually broadly follows topography [22]. Thus, by combining potentiometric surface maps, hydraulic cross-sections, pressure-depth diagrams, maps of hydraulic-head differences, and ecological site maps, or by incorporating numerical simulation, it becomes possible to estimate the discharge and recharge areas' spatial distribution in a region of interest [22,23].
At the European level, the Water Framework Directive (WFD; 2000/60/EC) emphasizes the role of parallel integrated quantitative and qualitative monitoring and evaluation to achieve a good general status for the main European water bodies, including aquifers. Consequently, the groundwater budgeting and vulnerability assessments of main aquifers should be analyzed using such integrated approaches [25,26]. Recently, automated shallow groundwater monitoring networks and large datasets have become available in most countries, and these have also been used to monitor groundwater depletion [9,27,28]. One agricultural area where SGW is the subject of both anthropogenic exploitation Black diamonds represent the SGW wells used in the current study, while those not used are shown in white. The digital elevation model (DEM) was constructed from SRTM data [45]. Hydrogeological type section of the Danube-Tisza Interfluve with the aquifer (c) based on Mádl-Szőnyi and Tóth, 2009 [37].

SGW Level Variation
In this study, 190 shallow groundwater (SGW) level time series (covering 1961 to 2010) of monitoring wells operated by the Lower Danube Valley Inspectorate (Figure 1a) were used. The average screening depth was 6.8 m (maximum: 20 m). Thus, every monitoring well can be said to Black diamonds represent the SGW wells used in the current study, while those not used are shown in white. The digital elevation model (DEM) was constructed from SRTM data [45]. Hydrogeological type section of the Danube-Tisza Interfluve with the aquifer (c) based on Mádl-Szőnyi and Tóth, 2009 [37].

SGW Level Variation
In this study, 190 shallow groundwater (SGW) level time series (covering 1961 to 2010) of monitoring wells operated by the Lower Danube Valley Inspectorate (Figure 1a) were used. The average screening depth was 6.8 m (maximum: 20 m). Thus, every monitoring well can be said to screen the uppermost, unconfined aquifer, which in most cases, consists of fine sand and silt (but no detailed information by well). The average density of the SGW well distribution was 2 wells per 100 km 2 ; their number and temporal sampling frequency, however, varied over the time interval studied (Figure 1b). The original data was monthly. Due to the annual focus of the present study, annual medians of SGW levels were calculated and only those years were considered in which at least 9 months were represented by data.

Indices of Meteorological Variability
Indices of meteorological variability were used as potential explanatory parameters in the study (see Section 3.3 for details). Annual data of average temperature ( • C), precipitation (mm y −1 ), potential evapotranspiration (mm y −1 ), sunshine (duration in hours-h), global radiation (J cm 2 ), the moisture index (-), the Palfai Drought Index (-), and the aridity index (-) were obtained from the CarpatClim database (1961 to 2010) with a grid resolution of 8 × 8 km [46,47]. A real average time series were computed for the area for the different parameters obtained from the CarpatClim database.

Long-Term "C umulative" Meteorological Indices
In addition to the gridded meteorological parameters (Section 2.2.2), cumulative meteorological variables were derived from precipitation and potential evapotranspiration, to better represent the elongated processes related to aridification prevailing in the area [32].
The Standardized Precipitation Index, (SPI) [48] is a well-known and widely used index for the determination of drought periods [49]. The calculation of the SPI is based on the long-term precipitation record of a given period (in this case 50 years). This long-term record is fitted to a probability distribution, which is, then, transformed into a normal distribution [48]. Values of the SPI can be expressed as extreme dry (SPI ≤ −2), severe dry (−2 < SPI ≤ −1.5), moderate dry (−1.5 < SPI ≤ −1.0), near normal (−1 < SPI < 1), moderate wet (1.0 ≤ SPI < 1.5), severe wet (1.5 ≤ SPI < 2.0), and extreme wet (2 ≤ SPI). The index was designed to quantify the precipitation deficit over multiple timescales [48]. In the present study, 12 and 24 month Standardized Precipitation Index were calculated (SPI−12 and SPI-24; no measurement unit) from monthly precipitation time series of the CarpatClim database with the precintcon R package [50,51].
The calculation of cumulative potential evapotranspiration departure (mm y −1 , CPED) and cumulative precipitation departure (mm y −1 , CPD) is a method used by climatologists and hydrologists to delineate precipitation trends [52] and is also sometimes used to describe groundwater recharge [53][54][55]. CPD and CPED time series were calculated as follows (Equation (1)) where x t is an element of the original time series at timestep t, y t is the element of the cumulative time series, andm the average of the original time series. Please note that SPI is calculated by using a moving window width that is equal to the number of months desired (in the present case 12 or 24 months), which smooths the high frequency fluctuations, thus, it is appropriate to study the quasi-and supra-decadal precipitation variance [48,56]. In the meanwhile, CPD and CPED cumulates the data of consecutive years rendering the dataset more applicable to study the low frequency variability and trends (full memory) practically without any high frequency signal [53].

Methodology
In order to determine the spatiotemporal change in the SGW table in the area, as well as to establish a picture of its driving parameters, the following steps were followed: After evaluating the spatial variability of (i), the SGW level minima and maxima (range), and (ii) the long-term temporal changes of water levels (trends), in step (iii) the groups of the SGW level time series were obtained using HCA. Finally, (iv) the linear relationship of the SGW level time series and their group means was evaluated to climatic parameters ( Figure 2). The set of methods were chosen based on their wide applicability in hydrogeology. Hierarchical cluster analysis combined with discriminant analysis has been used, for example, to group the SGW wells of a similar area in Austria [26] to try to delineate their common drivers [57], whereas variogram analysis was applied, for example, to determine and map the response of SGW time series to changes in recharge condition on the banks of the Danube [58]. However, this particular combination and order of application of these tools ( Figure 2) has not yet been used to delineate flow regime areas before.
Water 2020, 12, x FOR PEER REVIEW 6 of 24 using HCA. Finally, (iv) the linear relationship of the SGW level time series and their group means was evaluated to climatic parameters ( Figure 2). The set of methods were chosen based on their wide applicability in hydrogeology. Hierarchical cluster analysis combined with discriminant analysis has been used, for example, to group the SGW wells of a similar area in Austria [26] to try to delineate their common drivers [57], whereas variogram analysis was applied, for example, to determine and map the response of SGW time series to changes in recharge condition on the banks of the Danube [58]. However, this particular combination and order of application of these tools ( Figure 2) has not yet been used to delineate flow regime areas before.

Data Imputation and Aggregation
In the case of data gaps, the missing annual data were imputed using multivariate regression analysis computed between the neighboring SGW wells within a 50 km search radius if the number of complete cases (a requirement for hierarchical cluster analysis) between the SGW well pairs was larger than 30 (years). If their linear relationship was significant (p < 0.05) and r > 0.75, a maximum of two consecutive years was imputed; this occurred in 1.2% of the whole dataset. As a result, 91 SGW monitoring wells had data representing at least 80% of the time interval, while 73 of them represented the investigated 50 years continuously, while the other 98 SGW wells were not used in this study ( Figure S1).

SGW Levels Variation Indicators and Statistics
The first step in the analysis was to explore the descriptive statistics of the SGW levels, with a special focus on the spatial pattern of the range (maximum and minimum levels) of water level fluctuation within the investigated time interval (1961 to 2010, Figure 2). As a next step, linear regression models (Equation (2)) were used to calculate the linear trends of each of the SGW level time series using an ordinary least squares method (Equation (3)) where β0 and β1 are the coefficients of the linear regression model and ε is the random term.
where y is the t th estimated value of the trend line, and x is the t th measured SGW level at a given time, and b is the intercept b , the slope (the estimators of β and β ). The significant slopes (b1) were presented on interpolated maps derived using the methods discussed above. In simple linear

Data Imputation and Aggregation
In the case of data gaps, the missing annual data were imputed using multivariate regression analysis computed between the neighboring SGW wells within a 50 km search radius if the number of complete cases (a requirement for hierarchical cluster analysis) between the SGW well pairs was larger than 30 (years). If their linear relationship was significant (p < 0.05) and r > 0.75, a maximum of two consecutive years was imputed; this occurred in 1.2% of the whole dataset. As a result, 91 SGW monitoring wells had data representing at least 80% of the time interval, while 73 of them represented the investigated 50 years continuously, while the other 98 SGW wells were not used in this study ( Figure S1).

SGW Levels Variation Indicators and Statistics
The first step in the analysis was to explore the descriptive statistics of the SGW levels, with a special focus on the spatial pattern of the range (maximum and minimum levels) of water level fluctuation within the investigated time interval (1961 to 2010, Figure 2). As a next step, linear regression models (Equation (2)) were used to calculate the linear trends of each of the SGW level time series using an ordinary least squares method (Equation (3)) where β 0 and β 1 are the coefficients of the linear regression model and ε is the random term.
whereŷ t is the t th estimated value of the trend line, and x t is the t th measured SGW level at a given time, and b 0 is the intercept b 1 , the slope (the estimators of β 0 and β 1 ). The significant slopes (b 1 ) were presented on interpolated maps derived using the methods discussed above. In simple linear regression, the t-test and F-test are used to determine if the linear relation between variables is significant. Both tests lead to the same conclusion. The hypothesis of both tests evaluates if the value of β 1 is zero. If the null hypothesis is rejected, then it must be concluded that β 1 0. In the present case, the significance of the slopes was tested with the use of the F-test at α = 0.95 [59].

Grouping of the SGW Wells and Validation
Prior to spatial grouping, the SGW time series were Z-score transformed, i.e. overall average SGW level was subtracted from the distinct SGW level measurements in a given SGW well, and divided by the standard deviation of the water level in that particular well. The spatial grouping of these Z-score transformed SGW time series was assessed using hierarchical cluster analysis (HCA) following Ward's method [60] and squared Euclidean distance [61]. Then, the grouping was validated using linear discriminant analysis (LDA) [62] and the Dunn index [63]. The former was performed by classifying the nearest-mean type in the transformed space using the LDA function of the MASS package [59] in R [51]. The latter, however, is a metric for evaluating clustering algorithms. The higher the index value is, the more compact and well-separated the clusters are from each other.
HCA is one of the most frequently employed clustering algorithms in cases where the aim is to group objects on the basis of their measured parameters [64]. In the present case, agglomerative HCA, in which the SGW wells formed separate clusters, was applied, and at each step, the two closest clusters were iteratively merged into one cluster with their linkage distance measured as squared Euclidean distance, the better to separate the clusters. The number of clusters decreases by one with each step, until only one, single cluster remained. The results obtained were plotted on a dendrogram which was intersected at 60% of the maximum linkage distance. Depending on the intersection, different groupings may be obtained.
In the interest of validating the cluster results, the question of whether the maximal separation between the groups and the simultaneous minimization of the variability was achieved with the obtained grouping was raised. Fisher's linear discriminant analysis (LDA) was used in the transformations (i.e., linear combinations) of the original data [65] to answer the question, as was the case in numerous other studies in water science (e.g., [26,[66][67][68][69]). In the procedure, the original underlying class for each observation was derived from HCA, and the evaluation of the number of cases in which the LDA predicted the correctly classified memberships was expressed as a percentage. This, in turn, also functions as a measure of the degree of the separability of the classes [65].

Spatial Interpolation and Correlation with Meteorological Variables
The obtained groups of SGW wells were, then, characterized by correlating the groups' mean SGW levels with variables representing local climate variability, and also correlating the distinct SGW level time series belonging to a given group with the same set of external parameters (see Sections 2.2.2 and 2.2.3). The latter correlations were visualized on box-and-whiskers plots. Interpolated maps were derived from the correlation coefficients of the given SGW time series using kriging [70], in which the weights of the interpolation were obtained using variography [71]. In the present case, the empirical semivariograms indicate the changes occurring in the spatial autocorrelation of the data, statistics describing the SGW levels in space. Using this information, theoretical semivariograms were fitted to the empirical semivariogram of the observed data providing the numerical function to fit a model of spatial correlation to the observed phenomenon, which was used during the interpolation [71,72]. To obtain second-order stationarity [73] a quadratic spatial trend was removed from the direct input data to the variography.
In the interests of reproducibility, in the case of every interpolated map, the parameters of the underlying theoretical semivariograms are documented in the Supplementary Materials ( Figure S6). The most important characteristic of the variograms is their spatial range, which is the distance within which the samples have an influence on each other and beyond which they are uncorrelated [74,75], the type of the semivariograms, sill, nugget, number of bins, max lag distance, and the minimum number of pairs in a bin. For a detailed description of kriging with variography, the reader is referred to Webster, 2008 [71]. The semivariograms and the interpolation were executed in EOV (EPSG 23700), while for the sake of better visualization, the maps in the study are presented in Lat, Lon (EPSG 4326, see Figure 1).

Spatial Distribution of the Temporal Variability of the SGW Levels
In order to obtain a preliminary picture of the temporal variability of the SGW table of the area, the range (maximum to minimum), and standard deviation (st. dev.) of the SGW level of each well were determined for the years 1961 to 2010 (Table S1) and used to derive interpolated maps ( Figure 3). The greatest degree of variability in the SGW level, a range of 5 m, was observed in the central N-S "axis" of the area. Between the banks of the River Danube and the central part of the area, the range of SGW level variability was found to be quite small, approximately~1 m, while towards the River Tisza it increases to 1.5 m, indicating more stable SGW levels within the studied time interval. The st. dev. ( Figure S2) followed a similar linear pattern to the range (r = 0.98).
Water 2020, 12, x FOR PEER REVIEW 8 of 24 to Webster, 2008 [71]. The semivariograms and the interpolation were executed in EOV (EPSG 23700), while for the sake of better visualization, the maps in the study are presented in Lat, Lon (EPSG 4326, see Figure 1).

Spatial Distribution of the Temporal Variability of the SGW Levels
In order to obtain a preliminary picture of the temporal variability of the SGW table of the area, the range (maximum to minimum), and standard deviation (st. dev.) of the SGW level of each well were determined for the years 1961 to 2010 (Table S1) and used to derive interpolated maps ( Figure  3). The greatest degree of variability in the SGW level, a range of 5 m, was observed in the central N-S "axis" of the area. Between the banks of the River Danube and the central part of the area, the range of SGW level variability was found to be quite small, approximately ~1 m, while towards the River Tisza it increases to 1.5 m, indicating more stable SGW levels within the studied time interval. The st. dev. ( Figure S2) followed a similar linear pattern to the range (r = 0.98).

Determination and Characterization of the Spatial Groups of SGW Wells
After distinctive areas with different temporal behaviors in the SGW table were outlined (Section 3.1), the SGW levels were Z-score transformed and clustered (see in Section 2.3.3) to group the wells with similar behavior in terms of normalized SGW level variation. The results of HCA ( Figure S3a,b) indicated three sets of SGW wells with distinctively different behaviors (Figure 4a). This finding was also verified using LDA, indicating that the cases had been correctly classified with a 100% level of certainty ( Figure 4b); also, the Dunn index [63], gave the highest value when the SGW wells were split into three groups out of all the possible combinations ( Figure 4c).

Determination and Characterization of the Spatial Groups of SGW Wells
After distinctive areas with different temporal behaviors in the SGW table were outlined (Section 3.1), the SGW levels were Z-score transformed and clustered (see in Section 2.3.3) to group the wells with similar behavior in terms of normalized SGW level variation. The results of HCA ( Figure S3a,b) indicated three sets of SGW wells with distinctively different behaviors (Figure 4a). This finding was also verified using LDA, indicating that the cases had been correctly classified with a 100% level of certainty ( Figure 4b); also, the Dunn index [63], gave the highest value when the SGW wells were split into three groups out of all the possible combinations ( Figure 4c). To obtain a larger-scale picture, the trends of the SGW levels were investigated also at the group level, i.e., the slope of the areal group mean SGW level time series was also modeled. These also reflected a significant (p = 2.09 × 10 -9 ) drop (β1 = −0.042 m y −1 ) in water level in the wells of G1 (Table  1) from 1961 to 2010, with the water table reaching its maximum in 1966 and minimum in the mid-1990s (Figure 5c). An average decrease in water level also characterized G2, although this was one order of magnitude lower, and statistically insignificant (p = 0.0734), while there was a minor, but nonetheless significant (p = 3.85 × 10 −4 ) increase in SGW levels in G3 (Table 1). The three groups of the SGW wells were almost perfectly separated in space. The wells in Group 1 (G1, indicated by black triangles in Figure 4a) were clustered in the Central Ridge of the Danube-Tisza Interfluve at altitudes above 105 m aBSl, with the exception of eight wells scattered mostly in the eastern parts. This particular 105 m aBSl boundary falls within the lower quartile (q25) of the elevation of the SGW wells in Group 1 (G1), and the maximum in Group 2 (G2, indicated by white symbols), while no well was located at or above it in Group 3 (G3, indicated by grey symbols respectively in Figure 4a, Table 1). The wells in G2 and G3 are much closer (10 to 15 km) to one or the other of the rivers, while in fact, the SGW wells in G3 are effectively restricted to the banks of the Danube covering quite small range with respect to elevation (Table 1). Table 1. Descriptive statistics of the cluster groups of the SGW wells. The first row represents the statistics related to the elevation of the SGW wells in a particular group, while the other rows are related to the SGW levels of the SGW wells in a particular group. It is also clear from the descriptive statistics that the greatest degree of mean water level variability over the five decades in the elevated SGW wells occurred in G1, while the least was in G3 ( Table 1). The mean SGW levels of the wells in G1 varied by~3 m, and in G3 the figure was less than half of this (Table 1).
For each SGW well, the semicentennial trend of the SGW table was investigated by fitting a linear trend and evaluating its significance (Table S1). In the case of 69% of the SGW wells (n = 50), a significant decreasing trend (β 1 < 0, p < 0.05) was observed (Figure 4a), and 78% of these wells were clustered in G1 (n = 39, located in the Central Ridge), while 11 wells fell into G2, which was characterized by a less steep, yet still significant slope. There were 13 SGW wells (18% of the total) in which the slope of the water table change was insignificant within the studied 50 years; most of these (n = 11) belonged to G2, and only two wells belonged to G3 (Figure 4a and Table 1). In just a fraction of the SGW wells (13%), an increasing trend was observed. All of these belonged to G3 (n = 10) and could be delineated by the β 1 = 0.005 m y −1 contour line; i.e., a line representing the average increase in SGW level (Table S1), while the slopes in the wells of G2 were between −0.01 m y −1 < β 1 < 0.005 m y −1 , and wells of G1 were β 1 < −0.01 m y −1 (Figure 4a). The most significant drop in the SGW level occurred in the central parts of the study area, where the water level decrease even reached a value of β 1 = −0.11 m y −1 , equivalent to a 5.5 m drop over 50 years. This was also reflected in the greatest variability and range of water levels ( Figure 3 and Table 1). This negative trend, however, diminishes towards the rivers, while in the vicinity of the Danube, it even turns positive, indicating an increasing trend, with increases by 1 to 2 cm per year in SGW levels.
To obtain a larger-scale picture, the trends of the SGW levels were investigated also at the group level, i.e., the slope of the areal group mean SGW level time series was also modeled. These also reflected a significant (p = 2.09 × 10 −9 ) drop (β 1 = −0.042 m y −1 ) in water level in the wells of G1 (Table 1) from 1961 to 2010, with the water table reaching its maximum in 1966 and minimum in the mid-1990s (Figure 5c). An average decrease in water level also characterized G2, although this was one order of magnitude lower, and statistically insignificant (p = 0.0734), while there was a minor, but nonetheless significant (p = 3.85 × 10 −4 ) increase in SGW levels in G3 (Table 1).

Correlation Analysis with Meteorological Parameters
The observed differences in the long-term behavior of the SGW table in the different subregions (i.e., groups) called forth the investigation of possible external driving factors. In addition to the obvious anthropogenic drivers (e.g., water extraction and agriculture), the exploration of semicentennial common trends with meteorological variables (Sections 2.2.2 and 2.2.3) has been suggested as a focus of investigations [76,77]. Thus, first the areal mean SGW level time series of the groups was correlated with local meteorological parameters (Figure 5a), and this was followed by the SGW level time series of each well ( Figure S4 and Table S2).
Out of the available local meteorological parameters, the SPI, and the CPD indicated the three strongest meaningful relationships with the semicentennial change of the water levels in the groups with r > 0.7. In the case of G1, the CPD indicated the strongest linear relationship with the SGW level time series (r = 0.86), whereas for Groups 2 and 3 the SPI-24 gave the best result ( Figure 5a).
As a next step, the same investigative aim was pursued via the correlation of the SPI and CPD time series with the individual SGW wells' levels, looking at the spatial distribution of the correlations by group (Figure 5b). In the case of SPI-12 and SPI-24, quite similar results were obtained over the whole area. The wells in G2 and G3 indicate similar positive linear relationships with SPI, while the SGW level time series in the wells of G1 show a much weaker relationship with SPI in general (Figure 5b). On the contrary, CPD clearly determined the SGW variability of G1. This was because (i) SPI-12 and SPI-24 indicated a similar pattern, from which (ii) SPI-24 showed a stronger relationship with the SGW wells, and (iii) according to the literature, SPI-24 is very suitable for long term (i.e., semicentennial) investigations [78]; beyond this point, only SPI-24 was used, along with CPD.
The change in the cumulative precipitation departure (1961 to 2010) mirrored the change in the average of the SGW levels of the wells in G1 located in the central parts of the area (Figure 5c), while the linear relationships of the SGW time series with SPI-24 were weak. From 1965 to 1975, the value of CPD was positive, then, started to decrease intensively (β 1 = −0.039 mm y −1 , p = 1.85 × 10 −8 ) until 1995, as did the SGW levels. After 1995, CPD increased until 2010, while the SGW levels continued to decrease in a moderate way, with the lowest SGW levels observed in 2009.
Near the banks of the rivers Danube and Tisza (G2), the linear relationship of the SGW levels with the CDP was much weaker (Figure 5d), unlike with SPI-24. Near the rivers, the drop in the SGW level was not as severe as in G1; the lowest SGW level was in 1995, concurring with the semicentennial minimum of the CDP. The SGW table time series of G3 (Figure 5e), in the western part of the area, also had a strong connection with SPI-24, but because of the increasing trend of the SGW level in this particular subregion, the correlation was weaker than in G2 ( Figure 5).
It should be noted that SPI-24 displayed an insignificant positive trend during the investigated 50 years (β 1 = 2.32 × 10 -5 mm y −1 , p = 0.37), with a recurring significant periodic signal of 4 to 5 years (at p = 0.99; obtained with REDFIT . This particular 4 to 5 y periodic signal was commonly found to be significant against the red-noise background from an AR1 process (p = 0.8) in the areal average time series of all the groups but was most significant (p > 0.99) in G2. Periodic signals from the same band have been documented previously in the region, e.g., from precipitation [80], other hydroclimate proxy data archives [81], and SGW levels from the Danube-Tisza Interfluve [82] as well; and the common signal found reinforces the link between the behavior of the SGW wells near the rivers' banks and SPI-24. It must be said, however, that further investigations are required to determine its large-scale driving factors.
The spatial distribution of the correlation coefficients between the unique SGW level time series and SPI-24 ( Figure 6a) and CPD (Figure 6b) indicated a complementary spatial pattern ( Figure S6d,e). The linear relationship between SPI-24 and the SGW table are the highest near the banks of the two rivers (Figure 6a), especially on the side of the Danube. and severely wet conditions in years 1966, 2000, and 2006 (SPI-24 > 1.5), while between 1978 and 1996 drier periods were recorded (SPI-24 < 0). This interval concurred with the steepest section of the CPD. The lowest SPD value was observed in 1994, indicating severely dry conditions (SPI-24 < -1.5). This particular 4 to 5 y periodic signal was commonly found to be significant against the red-noise background from an AR1 process (p = 0.8) in the areal average time series of all the groups but was most significant (p > 0.99) in G2. Periodic signals from the same band have been documented previously in the region, e.g., from precipitation [80], other hydroclimate proxy data archives [81], and SGW levels from the Danube-Tisza Interfluve [82] as well; and the common signal found reinforces the link between the behavior of the SGW wells near the rivers' banks and SPI-24. It must be said, however, that further investigations are required to determine its large-scale driving factors.
The spatial distribution of the correlation coefficients between the unique SGW level time series and SPI-24 ( Figure 6a) and CPD (Figure 6b) indicated a complementary spatial pattern ( Figure S6d,e). The linear relationship between SPI-24 and the SGW table are the highest near the banks of the two rivers (Figure 6a), especially on the side of the Danube.  On the one hand, in those areas where a significant decrease in SGW table was seen (Section 3.2), the correlation with SPI decreases (Central Ridge (G1) above 105 m aBSl). On the other hand, those are the parts where the highest correlations with CPD (r > 0.8) are found (Figure 6b). The strength of the linear relationship between the SGW level time series and CPD decreases towards the banks of the rivers (G2). Reaching the River Tisza, the correlation coefficient declines to r = 0.6, while towards the Danube the pattern is more complicated. The decrease of the correlation is much bigger (r = −0.4). This coincides with the zone (G3) where the slope of SGW level change was positive (Figure 4a).

Comparison of the Delineated Subareas with the Previous Studies
The subareas delineated with HCA overlap to a great extent with those determined in previous studies, although those were found (i) on the basis of SGW variation from different time horizons, 1976 to 2003 (Figure 7a) [40] and 1970 to 2004 (Figure 7b) [38], and (ii) the SGW wells of the area grouped primarily on the basis of geomorphology [38,40]. Unfortunately, the numerical results of these studies are not available, thus, direct numerical comparisons were not possible, only empirical ones (Figure 7c).

Comparison of the Delineated Subareas with the Previous Studies
The subareas delineated with HCA overlap to a great extent with those determined in previous studies, although those were found (i) on the basis of SGW variation from different time horizons, 1976 to 2003 (Figure 7a) [40] and 1970 to 2004 (Figure 7b) [38], and (ii) the SGW wells of the area grouped primarily on the basis of geomorphology [38,40]. Unfortunately, the numerical results of these studies are not available, thus, direct numerical comparisons were not possible, only empirical ones (Figure 7c). From the assessment of the SGW table time series for a 30 y period, it becomes clear that Groups A and B2 [40] resemble the Central areas of the present study (G1), while Group B1 [40] covers the areas of G2 and G3, near the banks of the two rivers (Figure 7a,c). In another study spanning an interval of ~35 years and incorporating the SGW table time series of almost twice as many SGW wells as used in the present study, five subareas were determined, however their grouping was not validated numerically [38] (Figure 7b). These five subareas also show a significant overlap with the results of the present research (Figure 7c), specifically the fact that Group C can be considered to cover G3, while in most of its elements, Group E is equivalent to the SGW wells on the banks of the rivers (G2), and the other three groups can be considered to overlap with the Central Ridge of the Danube-Tisza Interfluve (G1) (Figure 7b,c).
In all of these previous studies, the groups found in the Central Ridge of the area have been characterized by a declining SGW table driven by geomorphology, the long-term lack of precipitation [38,39,83], water extraction, forestation, and the construction of channels [38]. Nevertheless, no hydrogeological assessment was conducted in relation to the groups and subareas obtained apart from the characterization of the Central Ridge, and the observation that the wells along the banks of From the assessment of the SGW table time series for a 30 y period, it becomes clear that Groups A and B2 [40] resemble the Central areas of the present study (G1), while Group B1 [40] covers the areas of G2 and G3, near the banks of the two rivers (Figure 7a,c). In another study spanning an interval of~35 years and incorporating the SGW table time series of almost twice as many SGW wells as used in the present study, five subareas were determined, however their grouping was not validated numerically [38] (Figure 7b). These five subareas also show a significant overlap with the results of the present research (Figure 7c), specifically the fact that Group C can be considered to cover G3, while in most of its elements, Group E is equivalent to the SGW wells on the banks of the rivers (G2), and the other three groups can be considered to overlap with the Central Ridge of the Danube-Tisza Interfluve (G1) (Figure 7b,c).
In all of these previous studies, the groups found in the Central Ridge of the area have been characterized by a declining SGW table driven by geomorphology, the long-term lack of precipitation [38,39,83], water extraction, forestation, and the construction of channels [38]. Nevertheless, no hydrogeological assessment was conducted in relation to the groups and subareas obtained apart from the characterization of the Central Ridge, and the observation that the wells along the banks of the River Danube (Group C-G3) are characteristically different from the others, something possibly caused by local discharge. The present study aims to fill this gap and place the obtained groups into a hydrogeological context.

The Relationship between the SGW Table, Surface Elevation, and Possible Meteorological Drivers
The border between the main subareas (G1) and (G2) is at 105 m aBSl, which is approximately in the middle of the basin elevation-wise. This concurs with the knowledge about gravitational flow systems and can be summarized as follows: (i) The differences in surface elevation are a primary control for subsurface waterflow, causing the groundwater table to follow topography [22,84,85]; and (ii) groundwater discharge is not restricted to the lowest point of valleys (e.g., stream channels), but also extends up to the entire lower half of a basin [86,87], in the present case, this being at 105 m aBSl. Thus, it is reasonable to observe that the behavior of the SGW table, in the area, mirrors the surface elevation of the SGW wells.
A positive linear relationship was observed between the range values of the SGW levels and topography (r = 0.76) in 58% of the cases (Figure 8a). This implies that as the elevation of the surface increases, so does the variability of the SGW level. A similarly "strong", but negative, linear relationship is present (r = −0.77) between the elevation of the SGW wells and the slope of SGW change in 60% of the cases (Figure 8b). This also implies that in elevated areas (Central Ridge, G1) the drop in SGW levels is bigger, unlike in less elevated ones. This is a key characteristic of gravitational flow systems, where the changes in the SGW table are more pronounced in the elevated and less so in the lower areas [22], as occurred here. These characteristics, however, did not account for the separation of G3 from G2. The meteorological parameters were applicable to answering that question.
Water 2020, 12, x FOR PEER REVIEW 16 of 24 example, compaction-driven or tectonic-driven flow [89]. The background of this explanation is explored in detail in Section 4.3. The four parameters' changes with elevation outline the differently behaving subareas to differing degrees; however, on the basis of the previous findings it is reasonable to suspect the following: • the Central Ridge of the study area (G1 with SGW wells located above 105 m aBSl) is fundamentally a recharge zone of a gravitational flow system (R); • although the area near the banks of the two rivers and parallel to them (G2) is characterized by In the Danube-Tisza Interfluve, a period with a significant lack of precipitation caused a continuous decline in the SGW up to the mid-1990s [38,39,88] (Figure 5c). This is illustrated on the SGW potentiometric maps from 1966 and 1995 ( Figure S5a,b). The two indices representing the lack of precipitation display characteristically different behaviors with respect to elevation (Figure 6a,b and Figure 8c,d). The SPI-24 correlations (Figure 8c) indicate a weaker negative linear connection with surface topography than the previous two parameters (r = 0.69), while CDP has the weakest (r = 0.58) and even an opposite (positive) linear relationship (Figure 8d). It should be noted, however, that the SGW wells assessed are screening the shallow, unconfined aquifers, thus, the change in precipitation amount will have an effect on the whole study area, regardless of the regime characteristic (recharge and discharge) of the subarea [22]. This is why the 4 to 5 y periodic signal is inherited from the climate parameters (SPI-24) by the SGW table, as documented previously in the area [82].
SPI-24 primarily represents interchanging wet and dry periods on a decadal scale [48] (Section 2.2.3). Thus, it was found to display a good linear relationship with the SGW table fluctuations of the lower-elevated subregions on the banks of the two rivers (G2 and G3); nevertheless, the~4 to 5 y fluctuation of the water table did not cause a significant overall change of the semicentennial timescale ( Figure 5d). This is a typical characteristic of the discharge zones, where the subsurface flow is able to counteract external drivers on a multi-annual scale (change in annual precipitation) to a certain extent [22]. This relationship was strongest in the lowest subareas and weakened towards the more elevated Central Ridge (G1) of the region (Figures 8c and 6a). This also occurred because as the elevation increased, subsurface recharge was unable to compensate for the lack of precipitation, drawing a boundary at the previously mentioned elevation of 105 m aBSl. This is where CPD becomes more dominant as an external driver, since it is primarily a low-frequency proxy for a long-term precipitation deficit [53] (Figure 5c), driving the SGW wells' table in the Central Ridge (G1). Therefore, it can be stated that the CPD index better describes the recharge areas, and the SPI-24 the discharge areas.
In addition, the behavior of the SGW wells in G3 does not seem to follow the long-term changes in precipitation (Figure 8d), and also displays a weak relationship with topography ( Figure 8d). The degree of the separation of G3 from the other two areas was even greater, to such an extent that would have been noticeable even without the application of HCA. This suggests the presence of an additional driving factor in the subarea of G3 on the banks of the River Danube, which could be, for example, compaction-driven or tectonic-driven flow [89]. The background of this explanation is explored in detail in Section 4.3.
The four parameters' changes with elevation outline the differently behaving subareas to differing degrees; however, on the basis of the previous findings it is reasonable to suspect the following: • the Central Ridge of the study area (G1 with SGW wells located above 105 m aBSl) is fundamentally a recharge zone of a gravitational flow system (R); • although the area near the banks of the two rivers and parallel to them (G2) is characterized by far less variability in SGW table than G1, it still has a connection to topography, so it seems likely to be the discharge zone of a gravitational flow system (D); • the SGW wells characterized by a slight increase in the water level located near the banks of the River Danube (G3, Figure 5e, Section 3.3) seem to belong to a discharge zone of another flow system (DO).

Embedding the Results in Hydrogeological Context
Previous studies have concerned themselves with the fact that in the study area, two characteristic flow regimes are present [90]. One is a gravity-driven flow system, which is controlled by the topography, the other one is an overpressured flow system developed by tectonic compression of the basement [37,90]. The former is a regionally unconfined flow system, characterized by normal fluid pressures recharged from precipitation in the Central Ridge of the area (Figure 9e). The increasing trends of SGW level (see Sections 3.2 and 3.3) experienced in the overpressured flow system (DO, G3) could originate in the alteration of subsurface hydraulic potential distribution caused by groundwater level decrease over great spatial and temporal scales. Groundwater recharge areas are more sensitive to changes in the amount of precipitation than discharge regions [22]. In the case of the latter, higher order groundwater flow systems could well buffer the effects of decreasing recharge, as reflected in Figure 9a, resulting in more stable groundwater levels. Because of the  (Figure 1a) provided by Mádl-Szőnyi and Tóth [37] (e). For further details on panel (e) see [37]. The color scale of the section lines (a-d) matches the original map color scale in Figure 3, Figure 4, Figure 6.
The other flow regime is an overpressured flow system developed by the tectonic compression of the basement. This regionally confined flow system is characterized by a dominantly vertical upward flow direction from the pre-Neogene units (Figure 9e) and the presence of a high concentration of total dissolved solids [37,90]. The hydraulic continuity of the two flow systems is facilitated by tectonic structures and highly permeable sandy lenses (Figures 1c and 9e) [37].
From the preceding results, it becomes clear that the elevated Central Ridge of the Danube-Tisza Interfluve acts as recharge area with dominantly downward flow direction (G1, R). The descending water of the gravity-driven flow system blocks the upwelling water of the overpressured flow system from reaching the surface (Figure 9e), as well as forcing the ascending deep waters to deflect towards the west and east. Along the discharge areas (G2, D) where flow is diverted upward by both flow systems, water with a large amount of total dissolved solids from the overpressured flow can reach the surface. However, the meteoric waters do not have sufficient momentum to push the deep waters westwards to the Danube, forming the smallest subarea (G3, DO) in the region. Consequently, the deep saline water is discharged before reaching the river, creating a saline tract in the area [91]. Discharge of the overpressured flow system is concentrated by tectonic structures and conductive formations (i.e., Gravel Aquifer, Figures 1c and 9e) in a distinct zone near the Danube [37,90].
These results can be compared to the ones of a hydrogeological study by Mádl-Szőnyi and Tóth [37] across a hydraulic section intersecting the study area ( Figure 9). According to the previous findings, the boundary between the gravity-driven and overpressured flow systems is determined by a characteristic difference in average annual change in SGW levels ( Figure 4). The SGW wells in the overpressured flow system had an average increasing annual SGW level higher than 0.005 m yr −1 . However, clarifying detailed groundwater flow patterns was beyond the scope of this study, in which currently used methods proved to be adequate to delineate areas with different hydraulic regimes (i.e., recharge and discharge zones), and also confirmed the border between the different flow systems (Figure 9), in agreement with the findings of Mádl-Szőnyi and Tóth [37].
Differences between the western and eastern discharge areas, which can be observed on the sections drawn from contour maps of the investigated parameters (Figure 9a-d), could well originate in two main factors. Due to the topographical conditions of the western discharge area, the steeper hydraulic gradient leads to a more intensive groundwater flow as compared with that in the eastern discharge area. Another aspect is the presence of a highly conductive formation close to the surface (see Lower Great Plain Aquifer and Gravel Aquifer in Figure 9e), which also tends to increase the dynamics of the flow system. As a consequence, changes are more intensive than in the case of the eastern discharge area, which is dominated by formations with a low hydraulic conductivity (Fluvial aquitard in Figure 9e). The increasing trends of SGW level (see Sections 3.2 and 3.3) experienced in the overpressured flow system (DO, G3) could originate in the alteration of subsurface hydraulic potential distribution caused by groundwater level decrease over great spatial and temporal scales. Groundwater recharge areas are more sensitive to changes in the amount of precipitation than discharge regions [22]. In the case of the latter, higher order groundwater flow systems could well buffer the effects of decreasing recharge, as reflected in Figure 9a, resulting in more stable groundwater levels. Because of the significant groundwater level decrease of the elevated Central Ridge of the Danube-Tisza Interfluve (dominated by gravity-driven flow system), and moderate changes in water level of the discharge regions (dominated by overpressured systems), it was possible for the topography of the groundwater level to experience considerable change. This could lead to the modification of flow dynamics (i.e., to decrease of flow intensity), and also to change of hierarchy of groundwater flow systems previously described, for example, in Havril et al. [92]. In the case of the study area, the reset of the upper fluid potential surface could lead to the expansion of overpressured flow system and, consequently, a slight increase in the above-lying SGW level.

Outlook on the Semicentennial Decrease of SGW
The 24-month Standardized Precipitation Index proved to be an efficient driver of the discharge areas, while the cumulative precipitation departure turned out to mirror the long-term change of the SGW level in the recharge areas accurately, since it represents the multiannual departure from the average and normal rainfall. This study draws attention to the severe change induced by precipitation loss in the water balance of a vulnerable recharge area. Since the SGW supplies the root zone with necessary moisture through capillary pressure, a decline in the SGW table breaks this connection, endangering the surface ecosystem. In the case of the Danube-Tisza Interfluve, the SGW table dropped an average of~2 m in the recharge zone in the last 50 years. This decrease could also be given as the reason for the decrease in area or complete drying-up of the shallow ponds in the central part of the study area, for example, Szappanszék Pond [42]. Furthermore, the landscape structure and habitat pattern of the fen and alkali vegetation has also started to change [41]. With persistent aridification (e.g., [32]), the system would become transient, and these unwanted changes could also take effect in the discharge areas near the two rivers, as has occurred in certain SGW wells in the gravitational discharge zone (decreasing trends, Figure 5d), or in the overpressured discharge zone (moderated increase in SGW level, Figure 5e).

Conclusions
The approximately 30-year-long aridification of the study area has manifested itself in a spatiotemporally variable way in its semicentennial shallow groundwater table. The documented difference in the behavior of the clusters of the SGW wells provides the foundation for the delineation of the different flow regimes (recharge and discharge areas) prevailing in the Danube-Tisza Interfluve combining the use of multivariate data analysis methods and geomathematics.
With the grouping of the SGW wells' semicentennial time series, three major subareas were delineated in the region and verified by the results of previous studies. The critical 105 m aBSl threshold found between the elevated Central Ridge and the areas near the banks of the two bordering rivers was in agreement with the findings in the literature and proved to be of assistance in the characterization of the subareas. A recharge zone was outlined in the central parts of the area lying above 105 m aBSl, with an average −0.042 m yr −1 decrease in water level over the past 50 years. The discharge zones (<105 m aBSl) parallel to the banks of the rivers bordering the area have also suffered a loss of water, although this is a whole order of magnitude smaller than in the recharge zone. In addition, the discharge zone of an overpressured flow system was also identified, in which the water level rose, despite the severe drop in precipitation.
The question of whether the areal mean SGW time series of the subareas, and the individual wells' water table time series display a similar relationship to external climatic variables (e.g., integrated precipitation, potential evapotranspiration, hours of sunshine and sunshine hours, etc.) was investigated. In the subareas characterized by discharge conditions near the banks of the rivers, the 24-month Standardized Precipitation Index was found to be a significant external factor, while in the more elevated recharge areas (Central Ridge), the cumulative precipitation departure proved to be a reliable significant driving factor.
The combined application of multivariate methods (HCA, linear regression, and LDA) and variography to the SGW data and additional environmental parameters have proven to be sufficient to delineate the borders between the different flow-regime areas (recharge and discharge zones, and an overpressured flow system), and to provide an overview of their relationship with meteorology and topography over a semicentennial scale, as well as determining their primary characteristics. The excess knowledge provided by this study could, in turn, serve as a benchmark for further preparation for the severe negative effects likely to come about as a result of climate change in the area [32]. Therefore, it is vital to explore those regions which are most affected by and vulnerable to climate change.