Regression Analysis and Risk Assessment of Groundwater Levels †

The present research work uses Residual Kriging to estimate the groundwater level of an unconfined alluvial aquifer, as well as the trend function. The ground surface elevation is used as auxiliary variable in the trend model. Indicator Kriging is applied to detect potential vulnerable locations. Classical variogram functions are applied to determine the spatial correlation of the measurements. The risk of hydraulic head to lie below a threshold value is significant, mainly at the South and North parts of the aquifer, where the lower values of groundwater level are estimated, indicating that these areas require intense monitoring to ensure the water resources availability.


Introduction
An accurate estimation of hydraulic head spatial variability is indispensable to develop an effective groundwater resources management plan, especially in arid regions.Kriging [1], a method based on linear interpolation of neighboring points combined with the mean square estimate error minimization, is one of the most exact interpolation schemes.The predicted value of an unsampled location is a weighted average of the neighboring data.The weights involved in the estimator are calculated using the distance between observation and estimation points and the spatial dependence function [2,3].
Different types of Kriging have been applied to the interpolation of hydraulic head [4][5][6][7][8].Ordinary Kriging (OK) is the most common variety of the homonymous interpolators family.The estimated values are based only on the sampled primary variable [9].Residual Kriging (RK) combines a trend function and the interpolated -by means of OK-residuals.The regression of the principal variable is depended on auxiliary information, usually related to topographic features.The deterministic trend model of groundwater level may incorporate space polynomials [10], elevation and topographic index [2,11], rainfall data [12].Indicator Kriging is used to estimate the probability of a variable to exceed or lie below a well-chosen threshold value at a specific point, mainly applied in risk analysis [13][14][15][16][17].
The present study applies RK to detect possible trend of the hydraulic head spatial variability of an alluvial aquifer, located at Northern Greece, considering in the estimator the ground surface elevation as secondary variable.Next, Indicator Kriging (IK) is used to predict the risk in terms of a set aquifer level limit (2.5 m above sea level), to ensure safe groundwater reserves.The two approaches provide useful information to develop an effective groundwater management strategy for the optimal manipulation of the area's ecosystem and agriculture.

Study Area
The proposed methodology is implemented in an unconfined aquifer that covers 210 km 2 , located at the Prefecture of Drama, Greece (Figure 1).The Water District of Eastern Macedonia, where the aquifer belongs in, is characterized by both continental and Mediterranean climate and the mean annual temperature varies from 14.5 to 16.0 °C.The annual rainfall ranges between 508 mm and 576 mm.
The observed data used to develop an accurate spatial model consist of 250 hydraulic head measurements.The available values vary between 0.9 m above sea level (masl) and 22 masl.Ground surface elevation, used as auxiliary information in the trend function, varies between 48 and 137 m.Higher values of elevation are met at North-East area of the aquifer.The dataset was provided by the Water Management Authority of Eastern Macedonia and Thrace and the Environmental and Spatial Planning Authority of the Prefecture of Drama.

Methodology
The aim of the proposed methodology is to derive estimates, using the values (in masl) at , 1, … , sampling points.An estimation point , lying inside the convex hall of the network, moves sequentially through all the nodes of the mapping grid.
To interpolate spatially the variability of groundwater level at the study area, RK is applied.A trend model that incorporates the ground surface elevation as auxiliary information in the drift term is developed, improving in this way the accuracy of the spatial estimation [11].The residuals are then interpolated with the help of OK.
In order to assess the risk of hydraulic head to lie below a threshold value, IK is implemented.The indicator transformed data points are obtained from a binary distribution and the conditional probability is calculated.
For each one of the two approaches, the experimental variogram is calculated and then fitted to classical theoretical models.To determine the suitable spatial model and evaluate the accuracy of the estimation, cross validation analysis is applied.The estimates are kept within the convex hull of the sampling points, as inside the quadrilateral are more accurate and precise than those outside.

Spatial Dependence
The variogram is commonly used in geostatistical analysis to express the spatial dependence of the monitoring points.The omnidirectional empirical (experimental) variogram is determined using the method of moments.The empirical variogram is fitted with isotropic classical models and the Matérn model, described in Table 1.

Table 1. Theoretical model functions:
is the variance, is the spatial separator vector, is the correlation length, is the coefficient, is the Hurst exponent, is the smoothness parameter, • is the Gamma function and • is the modified Bessel function of the second kind of order [18].
Model Equation

Residual Kriging
RK splits the random function into a smoothly varying and nonstationary trend and a random component, the residual data [9,19].Trend subcomponents can be determined based on physical principles.Herein, the deterministic trend of the hydraulic head is estimated by finding the relationship between the topographic elevation and the hydraulic head.The estimated values are expressed as a sum of the trend function and the interpolated by OK residuals ′ [20]: where : the deterministic regression coefficients, estimated using ordinary least squares (OLS) : the values of depended variables at the estimation point : the number of predictors : the kriging weights : the values of the sample data : the number of points included in the search neighborhood of the prediction point .
The variance of the estimation includes the trend prediction variance and the OK variance of the residuals: where : the vector of 1 1 predictors : the matrix of 1 1 predictors at the observations in the search neighborhood : the variogram matrix of the 1 1 residuals in the neighborhood : the Lagrange coefficient [21,22].

Indicator Kriging
IK is a useful geostatistical tool that helps to estimate the probability of a variable to exceed or lie below a well-chosen threshold value at a specific point, taking at the same time into account the uncertainty of the data [13].In the present research work, a limit that could cause problems of groundwater availability is set as a threshold value.IK estimates the probability of hydraulic head to lie below this value.The data are transformed into indicators using a binary distribution (e.g., 0 and 1).This transformation is given by the following equation: where is a binary variable, is the measured value and ′ is the cut-off (threshold) value.
The indicator variogram is estimated, using a suitable theoretical model.Kriging is applied using the usual equation to obtain predictions.The conditional probability at unsampled points based on the spatial correlation structure results in a map with values between 0 and 1, expressing the risk associated with the cutoff value [16,23].This probability map delineates suitable and unsuitable sites regarding the examined issue, while helps that way to take decisions to prevent and/or remediate a site compared to locations with reduced or no risk [24].

Results
As a first step, the data are decomposed into a linear trend plus a residual value.The experimental variogram of the residuals is then calculated and fitted to the theoretical models presented in Table 1.The optimal function that adapts the experimental omnidirectional variogram is the Exponential variogram.The estimated parameters are used to determine the spatial dependence model.To estimate the distribution of the hydraulic head data at unvisited locations and the associated uncertainty, RK is applied on a 100 × 100 grid.The corresponding contour maps are shown in Figure 2.
According to Figure 2a, the groundwater level changes from South-West towards North-East, following the ground surface elevation trend.The spatial variability of the groundwater level is established with higher accuracy around the sample locations, while at the North of the aquifer significant estimation error occurs (Figure 2b).To assess the risk associated with the groundwater level spatial distribution, IK is applied.As a threshold value is set the 25th, lower, percentile of the available data values, which is equal to 2.5 masl.The indicators are interpolated applying the best-fit Power-law variogram model.The probability of the groundwater level to lie below the cut-off value is increased towards the South and the North of the aquifer, where the lower hydraulic heads are met (Figure 3).

Conclusions
The Exponential variogram provides the best fit to the experimental one.The ground surface elevation, an auxiliary variable with physical notion, incorporates useful information about the hydraulic head spatial distribution, interpolated by means of RK.Due to the extremely low groundwater level values predicted at specific locations, IK is applied to identify the associated risk.The method is based on the best-fit Power-law model.The risk is significant mainly at the South and North parts of the aquifer, where the lower values of groundwater level are estimated, indicating that remedial actions should take place to avoid further decline of the aquifer.As a conclusion, RK delivers a quantitative estimation that shows how the groundwater level spatial distribution is formed, when other boundary conditions are not available.On the other hand, IK provides a measure of risk that the estimates calculated for unsampled locations by RK will be below the chosen critical threshold value.Overall, the application of the two aforementioned methods, accompanied by the appropriate spatial dependence function can provide useful and reliable results in order spatial variability and potential depletion risk of the aquifer level to be identified.

Figure 1 .
Figure 1.Map of the study area and its location in the Greek territory.

Figure 2 .
Figure 2. (a) Contour map of estimated hydraulic head values in the study area, using RK (red circles donate location of the 250 wells); (b) Contour map of estimated hydraulic head standard deviation using RK (blue circles denote location of 250 wells).

Figure 3 .
Figure 3. Map of estimated probabilities of the groundwater level of the aquifer to decline below the set 2.5 masl threshold value using IK.