Enhanced Characterization of the Krania–Elassona Structure and Functioning Allogenic Karst Aquifer in Central Greece

: The present study highlights the importance of geological, hydrogeological, and hydrogeochemical characterization of a karst aquifer in building a conceptual model of the system. The karst system of Krania–Elassona in central Greece was chosen for this application. Hydrogeological research included geological mapping and hydrogeological analysis. Additionally, hydrochemical analysis of water samples was performed in boreholes, rivers, and the system’s main spring. The Krania–Elassona aquifer consists of three horizons of marbles and is characterized by mature karstiﬁcation. The karst aquifer is characterized by allogenic recharge mainly from the River Deskatis that accounts for up to 92% of the total ﬂow. Groundwater and spring water are generally characterized as good quality and are suitable for irrigation and domestic use. The water type of the spring water is classiﬁed as Mg-HCO 3 . The application of a SARIMA (Seasonal Autoregressive Integrated Moving Average Model) model veriﬁed the conceptual model and successfully simulated spring discharge for a two-year period. The results of this study highlight the importance of basic hydrogeological research and the initial conceptualization of karst systems in reliably assessing groundwater vulnerability and modeling.


Introduction
Karst aquifers are a valuable source of fresh water worldwide [1] but are usually characterized by high vulnerability rendering them prone to external pollution [2]. Additionally, karst aquifers are susceptible to quantitative deterioration due to overexploitation [3]. However, they rarely coexist with human activities potentially harmful to groundwater and hence retain a good qualitative status. Nevertheless, increasing livestock activities and demands for agricultural land result in land use changes above karst hydrosystems and the abstraction of groundwater from these complex geological environments has increased along with their pollution risk.
In Greece, karst aquifers are of utmost importance as they sustain the water supply for agricultural, domestic, and industrial activities. Additionally, they significantly contribute to surface river runoff. However, in some cases surface runoff recharges the karst aquifers and contributes to the sustainability of these systems. The main karst hydrosystems are located in the western part of Greece where high levels of precipitation are observed and water demands are lower than in the eastern part of the The construction of SARIMA models requires stages of identification, estimation, and diagnostic checking [19,20,31]. The purpose of the identification stage is to determine the differencing required to produce stationarity and to estimate the order of both the seasonal and nonseasonal AR and MA operators (p, d, q, P, D, Q) for the stationary series. The identification is examined by the cumulative periodogram [C(fk)], the autocorrelation function (ACF), and the partial autocorrelation function (PACF) [19]. This information is then used to determine the general form of the univariate model to be fit. The estimation stage involves the estimation of the time series model parameters. This estimation is obtained by the residual sum of squares minimization using the Marquardt optimization algorithm for nonlinear least squares.  According to census data of the Hellenic Statistical Authority (2011), the permanent population of these settlements is about 14,500. Agriculture and livestock farming are the main economic activities occurring within the study area. Two rivers, Deskatis and Palaiomantanos, cross the north and south parts of the karst aquifer, respectively. The presence of the karst aquifer ensures the water supply and socioeconomic development of the study area. Moreover, the allogenic recharge of the karst aquifer is of special interest as such systems have rarely been investigated in Greece.

Hydrogeological Research
Determination of the hydrogeological regime of a karst aquifer constitutes the initial and possibly the most important stage to be performed before statistical analysis, hydrodynamic modeling, and vulnerability mapping. Hence, field work including geological mapping, measurement of tectonic features, hydrogeological mapping of karst features, groundwater water level measurements, groundwater sample collection, and pumping tests were carried out in the study area. The field work was supplemented with a detailed literature review and the evaluation of geoelectrical soundings and lithological profiles from previous studies [29], digitization of the geological and topographic maps in a GIS environment, creation of the site's digital elevation map, and analysis of the hydrochemical data using Aquachem and pumping tests.
Geological mapping was performed to update and verify the existing geological maps [30]. Geological mapping also helped identify the presence of marble horizons within the aquifer. Identification of tectonic features in the wider area is critical to determine the structure of the karst aquifer. Hence, faults, lineaments, and the stratification of the geological formations were identified to verify the existing data and/or include new features in the geological map.
Groundwater level measurements were taken from eleven wells, and the surface water flow of R. Deskatis was measured in three sites: the first located in the upper boundary of the aquifer, the second in the center of the system, and the third at the aquifer's outlet. Spring discharge was measured daily for four consecutive years (1989-1993) followed by monthly measurements until 2017. Hydrochemical analysis was performed in four groundwater samples, one from Kefalovriso spring and three from Deskatis River (Figure 1).
For the description of the seasonal SARIMA models (Seasonal Autoregressive Integrated Moving Average), the notations (p, d, q) (P, D, Q) S were used.
Where or with p = non-seasonal AR order, d = non-seasonal differencing, q = non-seasonal MA order, P = seasonal AR order, D = seasonal differencing, Q = seasonal MA order, and S = time span of repeating seasonal pattern.
Without differencing operations, the model could be written more formally as The non-seasonal components are: Note that on the left side of Equation (1) the seasonal and non-seasonal AR components multiply each other, and on the right side of Equation (1) the seasonal and non-seasonal MA components multiply each other.
Application of the SARIMA models requires stationarity of time series data obtained by different transformations. Logarithmic transformation is chosen to stabilize series variance and transform the usually skewed distribution into a normal distribution.
The construction of SARIMA models requires stages of identification, estimation, and diagnostic checking [19,20,31]. The purpose of the identification stage is to determine the differencing required to produce stationarity and to estimate the order of both the seasonal and nonseasonal AR and MA operators (p, d, q, P, D, Q) for the stationary series. The identification is examined by the cumulative periodogram [C(f k )], the autocorrelation function (ACF), and the partial autocorrelation function (PACF) [19]. This information is then used to determine the general form of the univariate model to be fit. The estimation stage involves the estimation of the time series model parameters. This estimation is obtained by the residual sum of squares minimization using the Marquardt optimization algorithm for nonlinear least squares.
The value of the parameters is associated standard errors, T-values, and P-values. Finally, diagnostic checking involves examination of the residuals fitted model which can, or cannot, prove the model's inadequacy and also provide information regarding model improvement. This determination can be achieved by using the identification stage tests and, furthermore, the Portmanteau test [32], the Akaike Information Criterion (AIC) [33] and the residual variance (σ e 2 ). For the log transformed series, the AIC criterions are given, respectively, as: where σ e 2 is the estimated residual variance, ln Z t is the sample mean of the logarithms of measured values, and n is the number of parameters. For a SARIMA model, the number of parameters n is given as:

Geological and Geomorphological Settings
The karst aquifer is crossed by the River Deskatis in the north and the River Palaiomantanos in the south. The hydrological basin of River Deskatis covers an area of 192 km 2 and crosses the karst aquifer for a distance of 19 km. The mean altitude is 759 m and the mean slope is 24.93%. The basin has a well-developed dendritic drainage network. The hydrological basin of River Palaiomantanos coincides with the karst aquifer over a relatively short zone (1000 m).
The study area belongs to the Pelagonian geotectonic zone, which consists of an underlying slightly transformed carbonate series known as the "Krania Carbonate Unity" [34,35] and the "Pelagonian tectonic cover". The latter tectonic cover comprises pro-alpine and alpine geological formations (gneiss, gneiss-schist, ortho-gneiss and marble layers) with plutonium boulders [36]. These formations are stone-paved in the Pelagonian zone and are covered by the Krania carbonate series. The "tectonic window" of Krania is the geological background of the Kam37vounion mountain range. It is distinguished by the appearance of three lithostromatographic horizons [30]: 1.
Lower horizon: The lower horizon occupies the eastern part of the tectonic window. It consists of medium-to thick-bedded, coarse to microcrystalline marbles, locally dolomitic, and is of light grey to white-grey color with a thickness of about 750 m.

2.
Mid horizon: The marbles of the intermediate horizon appear to be thin-bedded, locally foliated and rarely medium-bedded. The dark grey marbles of this horizon occupy the central part of the tectonic window. Their thickness reaches 700 m.

3.
Upper horizon: This horizon occupies the western part of the carbonate series and has a total thickness of about 400 m. It consists of medium-bedded and rarely thick-bedded marbles, locally dolomitic. They are mainly white in the lower parts of this horizon and grey at lesser depths. Very small percentages of lens-shaped mica schists occur in this section, which interfere with the crystalline marble of the top layer. The upper horizon is also characterized by sporadic mica schists. In many locations, such as NE of Loutro village, it is especially difficult to distinguish them from the expanded mica schists of the crystalline Pelagonian base.
The study site's tectonic cover is divided into two horizons. The lower horizon consists of ophiolite rocks (pyroxenes, serpentinites) in lenticular form, whereas the upper horizon is characteristic of Pelagonian cover and overlies the carbonate unity of Krania. Neogene deposits (marls, clays, and conglomerates) cover the tectonic drafts created by the tensile tectonics that took place from the Miocene onwards. The Quaternary formations comprise unlinked clay-sandy materials, pebbles, conglomerates, and gravel. The floodplain deposits of the Pleio-Pleistocene are sands and clays interfered with course materials, while the alluvial Holocene deposits are of river-basin origin and consist of sands, clays with loose pebbles, conglomerates, and gravel of various compositions.
The area is a tectonic basin and the main directions of the initial faults are NE to SW. The fault of Kefalovriso is located in the southeastern part of the study area. Hence, the fault adjacent to Kleisoura village indicates the NE boundary of the karstic aquifer [22]. In the broader area of Krania's tectonic window, only the tectonic structures created by the movements occurring during the alpine orogenetic phase of corrugations can be recognized [36]. Large anticline and convergent structures with curvature lengths of up to several kilometers can be distinguished. Their b-axes generally lie NW-SE and rarely ENE-WSW. According to the available geological maps of Elassona and Deskatis [30] and the field investigation, the following six groups of normal faults are present: (1) faults of NW-SE direction, (2) faults of NE-SW direction, (3) faults of NNE-SSW direction, (4) faults of ESE-WNW direction, (5) faults of NNW-SSE direction, and (6) faults of ENE-WSW direction. The second group (NE-SW) has played a significant role in the formation of the area's present geology, which is also due to the discontinuation of the Krania carbonates in the east.

Hydrogeological Setting
The karst aquifer of Krania-Elassona is characterized by mature karstification and consists of Mesozoic marbles. Three horizons of marbles occur within the aquifer. The upper horizon is highly karstified due to the tectonic stresses that occurred during the Pelagonian and that created adequate karstification conditions. The main karst forms are karren and dolines, the largest of which have diameters of up to 200 m and depths of 6 m. The distribution of the karst formations is basically controlled by the tectonic fractures with NW-SE and NE-SW direction. The high degree of marble karstification in the upper series was verified during drilling at depths of up to 150 m. The second horizon has similar karst formations. However, karst cavities are partially filled with red clay. In this layer, karstification reaches 300 m below the surface. The upper horizons recharge the lower horizon, while no springs are present in the first two horizons. In the lower horizon, surface karst formation rarely appears, and dolomite lines were detected during drilling. At greater depths of the lower horizon, the karstification degree is high. The aquifer is about 2000 m thick.
The central part of the karst aquifer is covered by Holocene sediments (alluvial and recent deposits), and the eastern part is covered by Neogene sediments. The aquifer has a mean elevation of 610 m and a perimeter of 62 km. The karst system is well defined by gneiss and schist outcrops located on its southern, northern, and western boundaries. The eastern boundary was delineated by using geoelectrical measurements [29]. The aquifer discharges from Kefalovriso spring. According to the pumping test, the transmissivity of the karst aquifer ranges between 1.1 × 10 −1 and 1.5 × 10 −1 m 2 /s, while the yield of the boreholes is up to 250 m 3 /h. Data from 1988 to 1993 showed that the piezometric head of the karst aquifer at this time varied from 255 to 257 m (Table 1). Groundwater flows mainly in a NW to SE direction ( Figure 2); however, on the eastern edge groundwater flows from north to south parallel to the aquifer's boundary ( Figure 3). The main characteristic of the system is its allogenic recharge from the River Deskatis. Figure 4 shows the monthly discharge of the river recorded from 1989 to 1993 in three sampling sites: (1) Deskatis bridge located before the aquifer's entrance, (2) Loutro site located 1.5 km from site 1 and within the karst aquifer, and (3) Valanida site located 12 km from site 1 and also within the karst aquifer. The three-dimensional graphic of Figure 4 presents the surface water losses that incur as they leave the crystalline rocks of the Pelagonian cover and continue to flow on the karstified marbles of the upper horizon of the Krania carbonate series. The total surface flow recorded in the five-year period 1989-1993 at the measuring point of Deskatis bridge was 20.5 × 10 6 m 3 . At the second sampling point of Loutro, the total amount of surface runoff recorded was 11.7 × 10 6 m 3 , while in the third site (Valanida), surface runoff did not exceed 1.8 × 10 6 m 3 . Percolation into the karst aquifer reaches 43% in the first 1.5 km of the river from its The total surface flow recorded in the five-year period 1989-1993 at the measuring point of Deskatis bridge was 20.5 × 10 6 m 3 . At the second sampling point of Loutro, the total amount of surface runoff recorded was 11.7 × 10 6 m 3 , while in the third site (Valanida), surface runoff did not exceed 1.8 × 10 6 m 3 . Percolation into the karst aquifer reaches 43% in the first 1.5 km of the river from its initial discharge, while at the 12 km point, up to 92% of the surface runoff percolates into the aquifer. It is worth mentioning that changes in these values are observed due to local contributions from tributaries and torrents of River Deskatis. In the south, River Palaiomantanos also crosses the karst aquifer for a distance of up to 300 m; however, it was not possible to obtain measurements in these sites. The aquifer discharges from the spring of Kefalovriso (Figure 1). The spring is a typical rising fault spring and is the only discharge point of the karst hydrosystem. Maximum discharge of the spring was recorded during the hydrological year 1990-1991 (3210 lit/sec), while the minimum discharge was recorded during the hydrological year 1989-1990 (259 lit/sec). The relatively low mean ratio between the minimum and maximum discharges recorded indicates a well-developed karst network within the aquifer. A representative hydrograph of Kefalovriso spring is presented in Figure 5. The hydrogeological analysis revealed the interactions between River Deskatis and the karst aquifer in terms of quantity. Nevertheless, rivers are usually recipients of nutrients originating from agricultural activities and these pollutants may influence aquifer quality.       The hydrochemical analysis was performed to determine the qualitative status of the groundwater, surface, and spring waters. The hydrochemical status of the system as well as the minimum, maximum, standard deviation, and mean values of the chemical analysis are shown in Table 2. The concentrations of calcium cations (Ca + ) observed in the three water types ranged from 30 to 68 mg/L, with 59 mg/L recorded at Kefalovriso spring. The concentration of Na + was significantly lower in groundwater samples than surface water samples with a mean value of 13 mg/L. The highest concentrations of SO − 4 were recorded in surface water samples as compared to groundwater samples. The highest concentration of NO − 3 (21 mg/L) was recorded in the surface water, while the concentration in Kefalovriso spring was 10 mg/L, highlighting the low impact of agricultural and livestock activities in this area. The concentrations of Mg 2+ ranged from 9.7 to 22.4 mg/L. The dominant water types are Mg-HCO 3 and Ca-Mg-HCO 3 . The Mg-HCO 3 water type can be attributed to Mg-bearing minerals found in dolomites, whereas the Ca-Mg-HCO 3 water type indicates groundwater of meteoric origin. The ion concentration plots of Piper, Durov, Schoeller, and Wilcox are presented in Figure 6. The ionic ratios of Ca 2+ + Mg 2+ /HCO − 3 , Mg 2+ /Ca 2+ were also calculated and are shown in Table 2. The Na + /Cl − ratio ranges from 1.0 to 2.9 indicating that waters are mainly meteoric in origin. According to the Wilcox plot, the samples are characterized by excellent-to-good quality and are thus suitable for irrigation.

SARIMA Model
SARIMA models were used in this research to simulate the discrete time series data of historical monthly discharge at Kefalovriso spring. The available historical monthly data cover twenty years (from January 1974 to December 1993) and are shown in Figure 7. The discharge measurements were obtained from the Land Reclamation Service (YEB) of Larissa and the Institute of Geology and Mineral Exploration (IGME). Of the available time series data, the period 1974-1993 was used to construct a suitable SARIMA model, and the years 1992-1993 were used as a forecasting period to check the model's reliability.

SARIMA Model
SARIMA models were used in this research to simulate the discrete time series data of historical monthly discharge at Kefalovriso spring. The available historical monthly data cover twenty years (from January 1974 to December 1993) and are shown in Figure 7. The discharge measurements were obtained from the Land Reclamation Service (YEB) of Larissa and the Institute of Geology and Mineral Exploration (IGME). Of the available time series data, the period 1974-1993 was used to construct a suitable SARIMA model, and the years 1992-1993 were used as a forecasting period to check the model's reliability. SARIMA models take into account both the seasonality, which was introduced by the length of seasons S, and the persistence of the time series. The stationarity of the time series data was obtained by logarithmic transformation. The Cumulative Periodogram [C(fk)], Autocorrelation Function (ACF), and Partial Autocorrelation Function (PACF) were computed for the logarithmic-transformed time series of monthly discharge (1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993) at Kefalovriso spring ( Figure 8). Figure 8 shows that the logarithmic-transformed time series of monthly discharge is not "white noise", the length of seasons S equals 12, and the model with structure (4,1,1) (1,1,1)12 is a candidate model. To allow for possible identification errors, a set of six models (Table 3) was considered. The number of parameters (n), the residual variance (σe 2 ), the Akaike Information Criterion (AIC), and the Portmanteau test (Qp) were applied to select the most suitable model (Table 3). Model (4,1,1) (1,1,1)12 had the lowest AIC and was selected to simulate the logarithmic-transformed time series of monthly discharge. As shown in Table 4, the model that satisfactorily fulfils most of the criteria is the (4,1,1) (1,1,1)12, which can be controlled [22].
The suitable seasonal stochastic model SARIMA (4,1,1) (1,1,1) 12 can also be used to predict the monthly discharge of Kefalovriso karst spring. The above model has the form: The suitable seasonal stochastic model SARIMA (4,1,1) (1,1,1)12 can also be used to predict the monthly discharge of Kefalovriso karst spring. The above model has the form: (  Therefore, SARIMA model (4,1,1) (1,1,1) 12 can be used to forecast monthly spring discharge for one or more time steps and was applied for the period 1992-1993. The observed and forecasted values of discharge at Kefalovriso spring were calculated using model (4,1,1) (1,1,1) 12 for the period 1992-1993 and are shown in Figure 10. The SARIMA model (4,1,1) (1,1,1) 12 predicted the spring's monthly discharge with sufficient accuracy. Therefore, SARIMA model (4,1,1) (1,1,1)12 can be used to forecast monthly spring discharge for one or more time steps and was applied for the period 1992-1993. The observed and forecasted values of discharge at Kefalovriso spring were calculated using model (4,1,1) (1,1,1)12 for the period 1992-1993 and are shown in Figure 10. The SARIMA model (4,1,1) (1,1,1)12 predicted the spring's monthly discharge with sufficient accuracy. Some basic statistics of the observed and forecasted discharge data at Kefalovriso spring for the period 1992-1993 are given in Table 5. The hypothesis that the mean and variance values of the forecasted data are not significantly different from those of the observed data can be accepted at the 5% significance level (Table 5). Thus, the results show that the forecasted data preserve the basic statistical properties of the observed data. Some basic statistics of the observed and forecasted discharge data at Kefalovriso spring for the period 1992-1993 are given in Table 5. The hypothesis that the mean and variance values of the forecasted data are not significantly different from those of the observed data can be accepted at the 5% significance level (Table 5). Thus, the results show that the forecasted data preserve the basic statistical properties of the observed data.
In this study, the allogenic recharge of the Elassona-Krania karst system was also quantified. The allogenic inflow contributes to the dissolution process introducing water that is undersaturated with respect to calcite [43,44]. Nevertheless, allogenic recharge should also be studied in respect to the diffuse recharge of karst aquifers [45,46] by using trace element analysis to distinguish allogenic runoff from karst waters. Future research in the karst system of Krania-Elassona should combine both hydrological modeling and hydrogeochemical analyses.    From the above analysis, it can be concluded that SARIMA models can be used successfully to predict short-term discharge values at Kefalovriso spring. Additionally, they can also generate sequences of synthetic monthly time series of unlimited duration. Therefore, their use can aid decision-making for the integrated protection and management of groundwater resources and lead to the appropriate programming of various management plans.

Research Challenges
The karst aquifer of Krania-Elassona is the main supply source of fresh water for the cities of Elassona and Krania, as well as the villages of Kefalovriso and Valanidia. Undeniably, the protection of the aquifer is crucial for the sustainability of the region. Thorough understanding of the geological and hydrogeological regime is the first step required to assess groundwater vulnerability and pollution risk of the Krania-Elassona system. According to the literature, various vulnerability methods have been applied to assess karst aquifers in Greece [37], including the COP methods applied to the Ziria system of Peloponnesus [38], the PaPRIKa method applied to the Palaiokastro-Chochlakies system on the island of Crete [39], and EPIK in Anthemountas basin [40]. However, in our opinion, statistical analysis and simulation of a karst system should precede groundwater vulnerability assessment. Statistical analysis has proven useful to determine the characteristics of a system [41]. Additionally, the large amount of high frequency data of the system could be used for rainfall-discharge analysis and modeling [42]. The determination of the conceptual model is the basis for a detailed model-based sustainable management plan for any karst system and its surrounding area. Kefalovriso spring discharges into the River Voulgaris, which is a tributary of the River Pineios. In this study, the allogenic recharge of the Elassona-Krania karst system was also quantified. The allogenic inflow contributes to the dissolution process introducing water that is undersaturated with respect to calcite [43,44]. Nevertheless, allogenic recharge should also be studied in respect to the diffuse recharge of karst aquifers [45,46] by using trace element analysis to distinguish allogenic runoff from karst waters.
Future research in the karst system of Krania-Elassona should combine both hydrological modeling and hydrogeochemical analyses.

Conclusions
The karst aquifer of Krania-Elassona located in central Greece is the main supply source of fresh water for the local community. The geotectonic evolution of the Hellenic region determined the tectonic window of Krania-Elassona. The main direction of the initial faults is NE-SW. The karst aquifer consists of three horizons of marbles and is characterized by mature karstification. Two rivers crossing the karst aquifer contribute to the recharge of the aquifer. The allogenic input from the River Deskatis accounts for up to 93% of the total runoff. The allogenic recharge from the Paliomantanos River is predicted to be much less than that of Deskatis. The karst aquifer discharges from Kefalovriso spring, which is a typical contact rising spring. The aquifer's groundwater quality is good, and the water type of the spring water is Mg-HCO 3 . Application of the SARIMA model verified the established conceptual model. Additionally, logarithms of monthly spring discharge time series can be simulated using SARIMA model. This model is suitable to simulate the Krania aquifer and can be utilised as a tool to predict monthly discharge values of Kefalovriso spring for, at least, a two-year period. Seasonal stochastic SARIMA models appear efficient to simulate both runoff and groundwater flow conditions of a karst system and can also be easily adapted to account for local conditions. Adapting a suitable stochastic model to specific karst groundwater flow conditions allows the user to obtain accurate short-term predictions that can aid rational groundwater resources exploitation and management planning.
The results of this study highlight the importance of basic hydrogeological research including geological mapping, detailed monitoring, and hydrochemical analyses, which constitute the basis for any model application. The successful application of SARIMA models enhanced system conceptualization and successful conceptualization constitutes the basis of reliable groundwater vulnerability assessment and modeling of any karst system.