On the optimal selection of interpolation methods for groundwater contouring: An example of propagation of uncertainty regarding inter-aquifer exchange

Accepted Manuscript

On the optimal selection of interpolation methods for groundwater contouring: An example of propagation of uncertainty regarding inter-aquifer exchange Marc Ohmer , Tanja Liesch , Nadine Goeppert , Nico Goldscheider PII: DOI: Reference:

S0309-1708(17)30614-0 10.1016/j.advwatres.2017.08.016 ADWR 2929

To appear in:

Advances in Water Resources

Received date: Accepted date:

13 June 2017 30 August 2017

Please cite this article as: Marc Ohmer , Tanja Liesch , Nadine Goeppert , Nico Goldscheider , On the optimal selection of interpolation methods for groundwater contouring: An example of propagation of uncertainty regarding inter-aquifer exchange, Advances in Water Resources (2017), doi: 10.1016/j.advwatres.2017.08.016

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights

AC

CE

PT

ED

M

AN US

CR IP T

9 interpolation methods for a Jurassic karst & a Quaternary alluvial aquifer are compared Calculated inter-aquifer exchange rates vary greatly depending on the chosen method Plausibility was additionally validated with eco-hydrogeological data (e.g. wetlands) Best results achieved with co-kriging incorporating additional data, e.g. topography

1

ACCEPTED MANUSCRIPT

On the optimal selection of interpolation methods for groundwater contouring: An example of propagation of uncertainty regarding inter-aquifer exchange 1*

1

1

Marc Ohmer , Tanja Liesch , Nadine Goeppert , Nico Goldscheider

1

* corresponding author 1

Institute of Applied Geoscience, Division Hydrogeology, Karlsruhe Institute of Technology (KIT), Kaiserstr. 12, 76131

Karlsruhe, Germany [email protected]

CR IP T

[email protected] [email protected] [email protected]

Abstract

The selection of the best possible method to interpolate a continuous groundwater surface from point data of groundwater

AN US

levels is a controversial issue. In the present study four deterministic and five geostatistical interpolation methods (global polynomial interpolation, local polynomial interpolation, inverse distance weighting, radial basis function, simple-, ordinary-, universal-, empirical Bayesian and co-Kriging) and six error statistics (ME, MAE, MAPE, RMSE, RMSSE, Pearson R) were examined for a Jurassic karst aquifer and a Quaternary alluvial aquifer. We investigated the possible propagation of uncertainty of the chosen interpolation method on the calculation of the estimated vertical groundwater exchange between the aquifers. Furthermore, we validated the results with eco-hydrogeological data including the comparison between

M

calculated groundwater depths and geographic locations of karst springs, wetlands and surface waters. These results show, that calculated inter-aquifer exchange rates based on different interpolations of groundwater potentials may vary greatly depending on the chosen interpolation method (factor 12.5). Therefore, the choice of an interpolation method

ED

should be made with care, taking different error measures as well as additional data for plausibility control into account.

AC

CE

PT

The most accurate results have been made with co-Kriging incorporating secondary data (e.g. topography, river levels).

2

ACCEPTED MANUSCRIPT Keywords: Groundwater contour map, Kriging, co-Kriging, Interpolation Method, inter-aquifer exchange, geostatistics

1. Introduction Reliable groundwater contour maps provide insight into manifold hydrogeological questions, e. g. determination of regional hydraulic gradients, flow directions, groundwater depths, flow velocities, recharge and discharge zones, hydraulic conductivities, aquifer susceptibility and catchment sizes to delineate protection areas. The comparison of contour maps of different points in time provides furthermore information about the temporal change and therefore the recharge and discharge of the area of interest. Differences in the hydraulic potential of two or more aquifers in the same area, separated by an aquitard, allow conclusions about possible vertical groundwater exchange.

CR IP T

With conventional methods, the groundwater level (GWL) can only be measured at distinct observation points, such as groundwater wells, springs and perennial surface water. Through the application of geostatistical and deterministic interpolation methods, the GWL can also be estimated between those measuring points. However, the number as well as the spatial and temporal distribution of the hydraulic head measurements are often, due to economic considerations, not sufficient to reliably represent the groundwater surface (Varouchakis and Hristopulos 2013; Delbari 2013). In hilly terrain, the estimation of the groundwater table (GWT) is often problematic because the data set is always sparse in relation to the topographic relief and measuring wells are almost exclusively located in valley region while the groundwater surface is usually a subdued replica of the ground surface elevation (Hoeksema et al. 1989).

M

AN US

Spatial interpolation methods and including geostatistics have been applied to various disciplines. A broad overview of comparative studies of interpolation methods in environmental science can be found in Li (2008). Zimmerman et al. (1999) compared the spatial interpolation accuracy of ordinary and universal Kriging and inversed distance weighting as well as the influence of surface types, sampling patterns, noise level and strength of small-scale spatial correlation on those methods by creating mathematical surfaces. They pointed out, that the Kriging methods outperformed the deterministic IDW methods over all levels and factors. During past decades, different types of univariate Kriging methods, e.g. ordinary Kriging (OK) and simple Kriging (SK) have been used to interpolate the GWT (e.g. Ahmadi and Sedghamiz (2007), Möhler et al. (2014), Sadat Noori et al. (2012), Guekie simo et al. (2016)) and were compared with one to several deterministic methods, e.g. inverse distance weighting (IDW) and radial basis function (RBF) (Sun et al. (2009), Hua et al. (2009), Chung and Rogers (2012), Yao et al. (2013), Varouchakis & Hristopulos (2013), Delbari (2013), Arslan (2014), Cooper et al. (2015) and Xiao et al. (2016).

PT

ED

In addition, some studies used multivariate kriging methods like universal Kriging (UK) and co-Kriging (CoK), to incorporate the influence of the topography on the interpolated GWT (Hoeksema et al. (1989), Ahmadi and Sedghamiz (2007), Sun et al. (2009), Chung and Rogers (2012), Sadat Noori et al. (2012), Delbari (2013), Varouchakis and Hristopulos (2013), Yao et al. (2013), Arslan (2014), Möhler et al. (2014), Cooper et al. (2015), Guekie simo et al. (2016), Xiao et al. (2016)). An overview of the methods and evaluation statistics used in those studies can be found in Tab. 1 and Tab. 2.

 

 

        

 

 

       









LR

CoCoK

DeK

                          

ANN

    

MC

                 

TSI

   

RK

                                                       

GWR

BK

CoOK

UK/KED

SK

OK

RBF

IDW

                                                                     

MLR

Ahmadi & Sedghamiz (2007) Arslan (2014) Chung & Rogers (2012) Cooper et al. (2015) Delbari (2013) Guekie et al. (2016) Hua et al. (2009) Möhler et al. (2014) Sadat et al. (2012) Sun et al. (2009) Tapoglou et al. (2014) Varouchakis & Hristopulos (2013) Xiao et al. (2016) Yao et al. (2013)

LPI

AC

Authors

GPI

CE

Table 1: Methods for interpolating groundwater contour lines applied in comparative studies. GPI: Global Polynomial Interpolation; LPI: Local Polynomial Interpolation; IDW: Inverse Distance Weighting; RBF: Radial Basis Functions; OK: Ordinary Kriging; SK: Simple Kriging; UK: Universal Kriging; KED: Kriging with External Drift; CoOK: Co-Kriging; BK: Empirical Bayesian Kriging; MLR: Multiple Linear Regression; GWR: Geographic Weighted Regression; TSI: Tension Spline Interpolation; MC: Minimum Curvature; DeK: Delaunay Triangulation; CoCoK: Collocated Co-Kriging; LR: Linear Regression; RK: Regression Kriging; ANN: Artificial Neural Networks; : Applied method; : Best method

    

             

   

   

 

 

 

 

   









3

ACCEPTED MANUSCRIPT

95 PPI

τ-Test

RMSSE

VSE

R²

                                                                                                                                                         

R

             

CR IP T

MSE

CVUD

RMSE

MAPE

MAE

                           

ME

Authors Ahmadi & Sedghamiz (2007) Arslan (2014) Chung & Rogers (2011) Cooper et al. (2014) Delbari (2013) Guekie et al (2015) Hua et al. (2009) Möhler et al. (2014) Sadat et al. (2012) Sun et al. (2009) Tapoglou et al. (2014) Varouchakis & Hristopulos (2012) Xiao et al. (2016) Yao et al. (2012)

OV

CV

Table 2: Evaluation statistics used in comparative studies. CV: Cross-validation; OV: Orthogonal-validation; ME: Mean error; MAE: Mean absolute error; MAPE: Mean absolute relative error; R: Correlation Coefficient; R²: Coefficient of determination; VSE: Variance of standardized error; CVUD: Cumulative Vertical Uncertainty standard deviation; MSE: Mean standardized error; RMSSE: Root mean square standardized error; τ-Test: Kendall’s rank correlation; 95PPI: 95 Percent Prediction Interval; : Applied statistics.



M

2. Methods

ED

 

Which interpolation technique provides the best results for the studied alluvial and karst aquifer? How do the different methods deal with the (i) spatially inhomogeneous distribution/patterns of the existing observation network? (ii) different surface types which change from a pronounced hilly topography to a flat riverine landscape? (iii) different hydraulic pressure conditions within the aquifer which fluctuates from unconfined to artesian? Which is the most suitable error statistics to compare the performance of the methods? How can the results be validated with additional eco-hydrogeological data? This includes e.g. comparison between calculated groundwater depth (GWD) and geographic locations of karst springs as well as wetlands and surface waters and finally comparison of calculated flow accumulation with the locations of receiving waters? What are the possible influences of the chosen methods on further computations, namely the calculation of the estimated vertical groundwater exchange between different aquifer systems?

PT

 

AN US

The majority of the previous studies compared only some of the methods and used selected error statistics for evaluation, and often did not investigate possible consequences of the chosen methods on further calculations or conclusions. In this study, almost all current deterministic and geostatistical interpolation methods, namely global polynomial interpolation, local polynomial interpolation, inverse distance weighting, radial basis function, simple-, ordinary-, universal-, empirical Bayesian and Co-Kriging, and error statistics, namely ME, MAE, MAPE, RMSE, R are systematically examined to answer the following research questions:

AC

CE

2.1 General overview For the purpose of this study, we employed a total of nine spatial interpolation methods on a Jurassic karst aquifer and a Quaternary alluvial aquifer. Of these were four deterministic methods (global and local polynomial interpolation, inverse distance weighting, and radial basis function) and six geostatistical methods (disjunctive-, ordinary-, simple-, universal-, empirical Bayesian-, and co-Kriging). For the co-Kriging method, additional surface-elevation data, river levels, and longterm groundwater levels of wells which were not measured at the reference date were included (see also 2.6.2). All model parameters were optimized by iterative cross-validation (CV). For clarity, the parameters are given in the supplementary. Spatial interpolation methods can be classified either as global or local methods. Global methods use all available data of the investigation area for the estimation and show a general trend, while local methods operate within a smaller area of variable size around the point to be estimated, and show a more local variation. Furthermore, the interpolation methods can be divided into exact and inexact methods. Exact methods generate a prediction equal to the observation at the sample point, while inexact (smooth) methods generate a prediction that usually differs from the observed values at the sample point (Li & Heap 2008), which may be appropriate for data that include considerable measurement errors compared to the total variations. Deterministic methods use closed-form mathematical formulas or the solution of a linear system of equations to estimate the value at a given location as a weighted sum of data values at surrounding locations, while geostatistical or stochastic methods include an assessment of a statistical spatial autocorrelation of the data by variography. In contrast to the geostatistical methods, which allow a probabilistic estimate of the interpolation quality, deterministic methods are not able to generate measures of uncertainty in addition to the estimations.

4

ACCEPTED MANUSCRIPT Almost all interpolation methods share the same formula to estimate the unknown value of ̂ at the point x0: ̂( )

∑

( )

(1)

Where Z(xi) is the measured value at the data point xi. The number of the existing data points is represented by n and represents the weight-function assigned to each data point (Oliver et al. 1996; Li & Heap 2008). 2.2 Deterministic Method

CR IP T

2.2.1 Global Polynomial Interpolation (GPI) GPI is a deterministic, global and inexact (smooth) trend surface analysis. It fits a smooth two-dimensional polynomial function of first, second or higher degree that represents the surface through a set number of data points (Cooper et al. 2015). There are a few decisions to make regarding the model parameters. The order of the function sets the shape of the interpolated surface. We achieved for both aquifers the best results with a first-order Polynomial which fits a flat plane through the dataset. The interpolated groundwater surface represents a gradual trend over the area of interest. The method is mainly recommended for groundwater surfaces that change slowly and gradually. Otherwise interpolated surfaces are highly susceptible to outliers especially in the edge region (Johnston 2004).

AN US

2.2.2 Local Polynomial Interpolation (LPI) Unlike GPI that adjusts a polynomial over the entire area, LPI adjusts several partially overlapping polynomials within a specified neighborhood. It is a moderately quick, inexact and local, but compared to GPI more flexible interpolator (Johnston 2004). The neighborhood shape, a minimum and a maximum number of points to be included, and the sector configuration can be specified. The surface value at the center of the neighborhood is estimated as the predicted value (Wang et al. 2014). It provides a prediction, prediction standard error and a condition number surface that are comparable to ordinary Kriging measurement errors. As with GPI, the order of the function determines the shape of the interpolated surface. As a constraint, the method is dependent on a regular distribution of the data points and normally distributed data values within the searching neighborhood (Johnston 2004). We used for both aquifers a first-order exponential Kernel function with 6 neighbors for the Quaternary and 29 for the Jurassic aquifer.

ED

M

2.2.3 Inverse Distance Weighting (IDW) IDW is an extensively used quick, deterministic, exact, and local interpolator. It estimates the value of a point by using a linear combination of the value at a sampled data point, weighted by an inverse function of the distance between these two points. The method is based on the assumption that observation points closer to the point of prediction are more similar to it as more distant points (Li 2008). The weights are computed according to the equation: ⁄ ∑

⁄

(2)

CE

PT

where di is the distance between the predicted point x0 and the measured point xi, n is the total number of measured points used in the interpolation and p is a power parameter which decides how the weight decreases as the distance increases (Xie et al. 2011). IDW cannot be used to predict values above or below the maximum and minimum measured values (Johnston 2004). With CV we determined a p of 4.88 for the alluvial and 10.26 for the karst aquifer. The shape of the searching neighborhood was chosen that all observation points will be used for the prediction.

AC

2.2.4 Radial Basis Function (RBF) RBF covers a large series of exact, moderately quick interpolators (e.g. thin-plate spline, spline with tension, multiquadric function, completely regularized spline) that are using a basic equation which is dependent on the distance between the predicted and the measured point (Aguilar et al. 2005). There is no assessment of prediction errors and they also do not allow investigating the auto-correlation of the data. The name comes from the fact that the function is radial-symmetric by definition. These functions are also used as basic functions of an approximation. RBF are a form of artificial neural networks (Chen et al. 1991) with a three-layer feedforward structure (input-, hidden-, output-layer). The input layer serves as an input distributor for the hidden layer. The hidden layer performs a fixed nonlinear transformation with no adjustable parameters and maps the input space into new space. The output layer implements a linear combiner on this new space and the only adjustable parameter is the weight of the linear combiner (Chen et al. 1991; Buhmann 2003).The approximants have the general form: ̂

( )

∑

(

)

(3)

5

ACCEPTED MANUSCRIPT where ̂ ( ) is approximated by a sum of n radial basis functions φ which have different centers ci and are weighted by the coefficients λi. We achieved the slightest RMSE with the spline with tension and completely regularized spline functions and have chosen the second one for both aquifers. 2.3 Geostatistical Methods 2.3.1 Kriging-Estimators Kriging is a very flexible interpolator that can be exact or smooth. It allows a variety of output surfaces including predictions, prediction standard errors, and probability (Johnston 2004). Each of the different Kriging methods is based on the following basic equation, which is a slight modification of equation (1): ̂( )

∑

( )

[ ( )

( )]

(4)

CR IP T

where λi is the Kriging weight derived from a covariance function or semivariogram, n is the total number of measured points, μ is a known stationary mean (trend component), assumed to be constant over the whole area of interest and calculated as the average of the data (Li & Heap 2011), and Z(xi) is the measured value at data point i. 2.3.2 Semivariance and Variogram Before the actual prediction, the spatial correlation of the data is assessed by variography. The semivariance γ of Z between the observation point xi and the prediction point x0 is defined as: (

)

( )

[ ( )

( )]

(5)

M

AN US

where h is the distance between the two points x0 and xi and γ(h) the semivariance. The plot of γ(h) vs. h is called empirical semivariogram and represents the spatial autocorrelation of the measured reference points. It quantifies the assumption that nearby data points tend to be more similar than more distant points. Some important features of the γ(h)-plot are the nugget, the sill/partial sill and the range. In theory, the semivariance at h=0 should be zero. The nugget is a positive value of γ(h) at h very close to 0. It represents variability at distances smaller than the typical sample spacing and includes measurement errors. The sill is the semivariance value at which the semivariogram levels off. The partial sill results from sill minus nugget. The range is the distance at which the sill is reached. Points that are further apart from each other than the range are considered spatially independent (Li & Heap 2011).

ED

There is a great variety of semivariogram models such as spherical-, exponential-, and Gaussian models, which in turn have significant influences on the prediction of the unknown values. The semivariogram model and the associated parameters nugget, sill and range are optimized in this study by using cross-validation with a focus on the estimation of the range parameter.

PT

2.3.3 Simple Kriging (SK) The estimation of SK is based on equations (4) and (5). It is assumed that the trend component is an exactly known constant over the whole area of interest and estimated by the mean value of measured data, μ(x0)= μ, so that: ̂ ( )

∑

( )

[ ( )

]

(6)

CE

We used for both aquifers a multiplicative skewing approximation method with a gamma distribution for the alluvial aquifer and a student t distribution for the karst aquifer.

AC

2.3.4 Ordinary Kriging (OK) OK is similar to SK with the difference, that μ is an unknown trend constant that has to be estimated. The most important consideration in OK is the assumption that the mean value remains constant over the whole are to be interpolated: ̂

( )

( ) ( )

∑

∑

( )

(7)

OK allows trend removal and data transformation. We interpolated the GWS with a smooth (OK Smooth) as well as a standard (OKStandard) neighborhood type. 2.3.5 Universal Kriging (UK) UK is also known as Kriging with a trend, Kriging with an external drift and regression Kriging (Hengl 2009). It is a multivariate extension of OK. Instead of a constant trend μ, it uses a linear or higher deterministic trend function μ(xi). A local trend function would be given by: ( )

(

)

(8)

6

ACCEPTED MANUSCRIPT As a trend function, we achieved the most satisfactory results with an exponential kernel function for the Quaternary and a Gaussian kernel function for the Jurassic aquifer. 2.3.6 Co-Kriging (CoK, CoOK) CoK uses information from one or more correlated secondary variables. The variables do not necessarily have to be measured at the same points but should be in the same range. CoK exists for all Kriging methods mentioned so far (where we used Co-ordinary Kriging (CoOK)). The main variable of interest Z1 and both autocorrelation for Z1 and crosscorrelation between Z1 and all other variables are used to improve the prediction. Cross-validation involves systematically removing single sample points from the data set, one after the other, and re-estimating their values by the used model (see also section 2.3). The prediction of Z1 cannot impair thereby because if there is no cross-correlation, it falls back to autocorrelation for Z1 (Rivoirard 1994; Johnston 2004).The secondary variables (support points) used for CoOK A-D are given in 2.6.2.

CR IP T

2.3.7 Empirical Bayesian Kriging (BK) BK provides a straightforward and robust interpolation method that automates the most difficult aspect of building a solid Kriging model by automatic calculation of parameters through a process of subsetting and simulation (Johnston 2004). It is using an intrinsic random function as the Kriging model. While other Kriging methods calculate the semivariogram from known data positions and use this single semivariogram model to make predictions for unknown positions, BK takes the error into account which is introduced by estimating the underlying semivariogram. 2.4 Validation Methods

M

AN US

2.4.1 Cross-validation (CV) The performance of a spatial interpolation method is affected by several factors, such as sampling density and distribution, and data variation (Li & Heap 2011). Therefore, the performance must be carefully evaluated in each case. The crossvalidation method is a statistical method to assess the accuracy of the interpolation. In cross-validation, each measured point is sequentially omitted, and the value is predicted by using the rest of the data. The difference between each measured and the respective predicted value is the error. Cross-validation can also be used to select the best possible modeling settings for the respective method (e.g. search radius, power option, kernel parameter). Based on the results of the cross-validation, following evaluation statistics or error measures were used to compare the accuracy of the different interpolation methods, where mi is the measured value and pi the predicted value at position i.

(

)

PT

∑

ED

2.4.2 Mean Error (ME) The ME is the average (arithmetic mean) of the errors. It indicates the average direction of the errors. An overestimation is indicated by positive bias, an underestimation is indicated by negative bias. ME is only conditionally suitable as an indicator of accuracy because negative and positive estimated counteract each other and the resultant ME tend to be lower than the actual error (Li & Heap 2011): (9)

CE

2.4.3 Mean Absolute Error (MAE) The MAE is the arithmetic mean of the absolute error values. It indicates the magnitude of the error which shows the accuracy of the method: ∑

(10)

AC

2.4.4 Mean Square Error (MSE) The MSE is given by the squares of the errors. It measures the magnitude of the error (accuracy), weighted on the squares of the errors, therefore it is very sensitive to outliers as it places a lot of weight on larger errors: ∑

(

)

(11)

2.4.5 Root Mean Square Error (RMSE) The RMSE is the square root of MSE. It has similar properties of MSE but has the advantage that it possesses the same unit as the required value (e.g. meter) [ ∑

(

) ]

⁄

(12)

2.4.6 Root Mean Square Standardized Error (RMSSE) The RMSSE should be close to 1. An RMSSE greater than 1 means a general underestimation in the variability of the predictions, an RMMSE smaller than 1 means a general overestimation in the variability.

7

ACCEPTED MANUSCRIPT √ ∑

[

]

(13)

2.4.7 Mean Absolute Percentage Error (MAPE) The MAPE express the accuracy as a percentage of the error: ∑

|

|

(14)

2.4.8 Pearson R and R² The Pearson R measures the linear correlation between the predicted and the measured rescaled covariance. Generally, it is in the range of -1 to 1. ∑

∑

Pearson

(15)

(∑

)

]√[

(∑

)

]

CR IP T

√[

∑

Often, its squared value R², known as the coefficient of determination, is given instead of R. Though R or R² are often used as a measure of the performance of the spatial interpolation methods, as discussed by Li and Heap (2008), it can be misleading and should be used with care.

ED

M

AN US

2.5 Study Area The study area is located in southern Germany, on the border of the federal states of Baden-Wurttemberg and Bavaria (Fig. 1). It comprises one of the most important sources of freshwater in Germany, which holds about one billion cubic meters and provides over 3 million people with high-quality drinking water (Flinspach et al. 1997). The region is hydrogeologically characterized by a complex interaction between the karst aquifer of the up to 550 m thick Upper Jurassic-Limestones and the alluvial aquifer situated within the Danube valley with its numerous ecologically important wetlands and the Danube River itself. The riverine landscape of the Donauried covers nearly half of the investigation area which is 2799 km². North of the Danube River, the rocks of the Upper Jurassic make up the high plateau of the Swabian Alb, which submerges below the Tertiary sediments of the pre-alpine Molasse with an average of 1-2° in SSE-direction. The low-permeable heterogeneous layers of the wedge-shaped Molasse might act as a hydraulic barrier and prevent a significant exchange between the karst and the alluvial aquifer. The north-eastern part of the study area is covered by large parts of the Miocene meteorite impact crater, the Noerdlinger Ries. The impact destroyed extensive areas of exposed rocks down to the crystalline bedrocks though, and only small-scale, economically unimportant groundwater reservoirs prevail (Winkler 1972). Therefore, it has been excluded for the interpolation of the groundwater levels.

CE

PT

2.5.1 Upper Jurassic karst aquifer As a result of its thickness, the Upper Jurassic consists of several individual karst aquifers. The groundwater flowing from the north to the Donauried belongs to the zone of deep karst. This deep karst zone is further divided (Villinger 1977) into an open deep karst zone and a zone of covered karst below the Molasse. In the transitional area between these two zones, the covering layers are either only partially present or thin so that the water level is unconfined in this northern edge region of the Molasse. The main karst aquifer is situated in the 200 m – 450 thick sequence above the Lacunosa-Marls (ki 1) starting with the advanced karstified Upper Kimmeridge Limestones (Felsen- and Bankkalke, ki 2-3) to the Hangende Bankkalk Formation (ti H). Regional layers of the impermeable Cement-Marl can lead to local areas with perched groundwater.

AC

2.5.2 Alpine Molasse The Paleogene/Neogene Molasse has intensive interbedded strata of clay to sandy layers. Accordingly, a distinction must be drawn between strongly anisotropic horizontal and vertical hydraulic permeabilities. The horizontal permeabilities of the -10 -4 individual layers are in the range of 5 x 10 m/s for silty clays and 4 x 10 m/s for silty sands. A large-scale effective vertical hydraulic conductivity should be rather at the lower end of this range. 2.5.3 Quaternary alluvial aquifer The relatively homogeneous fluviatile Quaternary gravels with their high to very high permeabilities represent very precious groundwater reservoirs. They cover the Molasse sediments completely within the Danube valley. They are found locally at the northern edge of the Molasse directly on top of the Jurassic limestones where a large-scale exchange of groundwater -4 -2 between the two aquifers can be assumed. Hydraulic conductivities between 5‧10 to 1‧10 m/s were determined in pumping tests (Bierer 1987).

8

ACCEPTED MANUSCRIPT 2.6 Available data 2.6.1 Water-level measurements Although a large number of measuring wells are available, there is a large discrepancy between the period and the interval of the measurements of the individual wells. The groundwater level measurements of March 11, 2013, were used as the basis for this study since this date ensured the best compromise between quantity and qualitative spatial coverage of the observation wells. Thereby, data from a total of 104 observation wells within the Quaternary and 77 within the Upper Jurassic (65 unconfined, 12 confined to artesian) were included for all methods.

CoOK A: monitoring wells with long-term values (LTV) (not measured at the reference day), CoOK B: digital elevation model data (DEM), CoOK C: DEM, river levels, CoOK D: LTV, DEM, river levels.

AN US

   

CR IP T

2.6.2 Support points Besides the reference date measurements, additional variables were used for co-Kriging interpolation to improve the prediction in areas with a low spatial coverage of data points. For areas with no available wells with measurements at March 11, 2013, long term values of wells measured in other years have been used as long they have a coherent measuring interval of at least 5 years and no observable trend. This results in additional data from 40 observation wells in the karst aquifer and 57 observation wells in the alluvial aquifer. A DEM based on SRTM1-data with a resolution of 30m x 30 m (1 arc second) was used to define the topography (USGS 2014).The surface elevation within the study area varies from 376 to 689 m asl. In addition, 16 river levels (daily average 11. March 2013, measuring stations: 13 Quaternary, 3 Jurassic) were used, where a direct contact of groundwater and surface water was expected. The following support points were used for the respective interpolation methods:

M

2.7 Data Processing The groundwater level data were evaluated using ArcGIS (10.3), Spatial Analyst and Geostatistical Analyst tools. The model parameters for each interpolation method, e.g. nugget, partial sill and others, were optimized using cross-validation with the focus on the estimation of the range parameter. A complete overview of all used model parameters is available in the supplementary material.

ED

2.8 Calculation of Groundwater exchange The large-scale hydraulic impact of the Molasse on the water exchange between the two aquifers has been little explored so far. The karst water level underlying the Molasse is in any part confined to artesian. A higher groundwater potential of the karst groundwater compared to the groundwater in the porous aquifer suggests areas with a possible infiltration upwards into the porous aquifer, whereas a lower potential a possible infiltration downwards into the karst aquifer. Provided a certain permeability is present on the spot

(13)

AC

CE

PT

The exchange of groundwater between the two aquifers in the area of the entire contact zone was quantified according to Darcy’s law by using the following equation:

9

ACCEPTED MANUSCRIPT -7

CR IP T

where K is the mean hydraulic conductivity for the leaky confining Molasse with 10 m/s. ∆h the potential difference between the two aquifers, A the size of the contact area and M the thickness of the Molasse which was determined by interpolation of data from 126 borehole cores.

AC

CE

PT

ED

M

AN US

Figure 1: Geological map of the investigation area

10

ACCEPTED MANUSCRIPT 3. Results and discussions 3.1 Cross-validation results Tab 3. shows the cross-validation results for all methods used in this study including error ranking. By combining the results of all validation methods in an average error ranking the order of suitability of the methods are for: the Jurassic aquifer: CoOK D > LPI > BK > CoOK A > CoOK C > CoOK B > OKSmooth > UK > RBF > OKStandard > GPI > IDW > SK. and the Quaternary aquifer: CoOK D > CoOK C > CoOK A > OKStandard > OKSmooth > BK > LPI > CoOK B > IDW > RBF > SK > UK > GPI.

CR IP T

Therefore, CoOK D, which uses additionally long-term values of monitoring wells not measured at the reference date as well DEM data and river levels, provides the best results for both aquifers. BK also provides good results for both aquifers. LPI achieves good results for karst while underperformed for the alluvial aquifer. The negative bias of the ME of LPI, CoOK A and CoOK D for the Jurassic Aquifer and LPI, SK, UK and CoOK D for the Quaternary aquifer indicates an underestimation while the other methods tend to overestimate the GWL. The CoOK methods achieved the best results for both aquifers. Within the co-Kriging methods, CoOK B (additional DEM values only) has the largest error parameters.

AN US

SK generates the greatest errors for the Jurassic and GPI for the Quaternary. The error estimations for the most frequently used methods IDW and OK show a substandard position in the ranking. The ranking of the methods shows no major deviations due to the use of different error statistics. Only Pearson R shows deviations for some methods which can lead to a different interpretation as long the method is compared with only one or a few other methods.

BK

CoOK C

CoOK D

BK

4 6 6 5 7 6

5 4 4 6 6 5

1 2 2 3 4 2.75

2 3 3 2 5 3.25

UK

CoOK A

CoOK B

CoOK C

CoOK D

BK

0.45 1.50 7.29 2.70 0.74 0.34 0.996

SK

-0.09 1.02 2.28 1.51 1.25 0.23 0.999

OKSTD

0.01 1.09 3.15 1.77 0.68 0.25 0.998

OKSMTH

0.34 1.53 9.08 3.01 0.68 0.35 0.995

LPI

0.26 1.41 4.62 2.15 1.28 0.32 0.997

RBF

BK

CoOK D

CoOK B

8 5 5 1 8 4.75

CoOK D

CoOK C

0.85 6.38 160.0 12.65 0.97 1.37 0.934

CoOK C

CoOK B

-0.47 6.38 149.8 12.24 1.49 1.38 0.938

CoOK B

CoOK A

0.26 6.65 167.8 12.95 0.73 1.44 0.932

GPI

0.25 -0.14 -0.16 1.37 2.33 3.12 4.75 16.10 18.79 2.18 4.01 4.33 0.95 0.25 3.57 0.31 0.54 0.70 0.997 0.992 0.990 Interpolation Method

0.74 6.62 188.7 13.74 0.60 1.42 0.922

CoOK A

UK UK

0.27 1.46 4.86 2.20 0.98 0.33 0.997

SK

-0.24 1.89 7.31 2.70 1.06 0.43 0.996

SK

OKSTD

0.27 2.19 15.55 3.94 0.50 0.991

OKSTD

OKSMTH

OKSMTH

LPI

0.04 3.26 18.71 4.33 0.74 0.990

-1.10 7.07 187.1 13.68 0.59 1.29 0.922

CoOK A

UK

LPI

ED

RBF

7 13 10 9 13 7 9 13 8 8 13 10 10 13 9 9 13 8.5 Interpolation Method

GPI

3 8 7 4 11 7.5

IDW

CE

1.06 2.57 1.34 6.93 14.06 7.81 198.3 492.2 196.6 14.08 22.18 14.02 2.13 0.66 1.52 1.47 3.12 1.67 0.918 0.867 0.919 Interpolation Method

IDW MAE MSE RMSE MAPE Pearson R Average

6 1 1 7 3 3

0.3 2.04 13.59 3.69 0.47 0.993

AC

Error Ranking

Validation Method Error Ranking

ME MAE MSE RMSE RMSSE MAPE Pearson R

11 11 11 11 2 8.75

SK

OKSMTH

RBF

12 12 12 12 1 9.25

PT

9 10 10 9 12 10.25

Validation Method

MAE MSE RMSE MAPE Pearson R Average Quaternary

OKSTD

LPI

1.14 6.52 196.6 14.02 2.84 1.39 0.918

M

RBF

-0.95 6.72 97.0 9.85 1.09 1.46 0.961

1.55 7.51 206.7 14.38 1.61 0.914

GPI

0.86 8.89 243.0 15.59 1.91 0.985

IDW

0.21 12.48 291.0 17.06 2.6 0.999

ME MAE MSE RMSE RMSSE MAPE Pearson R

IDW

GPI

Table 3: Cross-correlation. IDW: Inverse Distance Weighting; GPI: Global Polynomial Interpolation; RBF: Radial Basis Functions; LPI: Local Polynomial Interpolation; OK: Ordinary Kriging; SK: Simple Kriging; UK: Universal Kriging; CoOK -Co-Ordinary Kriging; BK: Empirical Bayesian Kriging; ME: Mean error; MAE: Mean absolute error; MSE: Mean Standard Error; RMSE: Root Mean Standard Error; RMSSE: Root mean square standardized error; MAPE: Mean absolute relative error; Pearson R: Pearson Correlation Coefficient. Error Ranking: 1 = “best method…13 = “worst method” Jurassic Interpolation Method

9 9 9 9 9 9

13 12 12 13 13 12.6

10 10 10 10 11 10.2

8 7 7 8 7 7.4

5 5 5 5 5 5

3 4 4 3 4 3.6

11 11 11 11 10 10.8

12 13 13 12 12 12.4

4 3 3 4 3 3.4

7 8 8 7 8 7.6

2 2 2 2 2 2

1 1 1 1 1 1

6 6 6 6 6 6

11

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Fig. 2 shows an error histogram grouped in 1 m (Quaternary, yellow) and 2.5 m groups (Jurassic, blue) as well as a scatter-plot of observed and predicted values. For the Jurassic aquifer, all methods tend to overestimate the low positioned groundwater levels in within the deep karst and to underestimate the high positioned water levels. This effect is most evident in SK and would result in a much lower gradient than in reality. Underestimated levels are primarily located in the area of the Swabian Alb. Those miscalculations could be the result of the variable topography or possible local areas with perched groundwater. Overestimated levels are located in the deep karst as well as in the western part of the investigation

Figure 2: Histograms for mean error (ME) and scatter plots of measured vs. predicted groundwater levels. CoOK: Co-Ordinary Kriging; OK: Ordinary Kriging; IDW: Inverse Distance Weighting; GPI: Global Polynomial Interpolation; SK: Simple Kriging; RBF: Radial Basis Functions; LPI: Local Polynomial Interpolation; UK: Universal Kriging; BK: Empirical Bayesian Kriging; [a]: best (CoOK D), worst (GPI) and most popular methods (IDW, OK) for Quaternary ; [b]: best (CoOK D), worst (SK) and most popular methods (IDW, OK) for Jurassic; [c]: other methods for Quaternary and [d] for Jurassic.grey dashed line: linear measured/predicted value; blue dashed line: R² predicted values.

area. Both are areas with a low density of the observation network. 3.2 Groundwater contour maps Though the comparison of the different interpolation results by error statistics shows differences, the overall differences for 2 2 example for the MAE or R are rather small, and, e.g. an R > 0.9 for nearly all methods suggests that all of them produce a good interpolation result. Therefore, additional methods for plausibility control are required, e.g. the inspection of the resulting groundwater contour maps. On the basis of to the pronounced anisotropy and heterogeneity of the karst aquifer due to the local alternation of hydraulic permeabilities, the interpolated contour maps can reproduce only a general large-scaled image of the groundwater surface.

12

ACCEPTED MANUSCRIPT 2

This can also be observed in the higher errors and lower R values compared to the Quaternary aquifer. In all contour maps, a basic flow direction of the karst water flowing from the Swabian Alb to southeast direction is shown. With the submergence of the Jurassic layers under the Molasse, the gradient decreases sharply in all interpolated maps especially in the western areas where many karst springs are situated at the boundary to the alluvial. Most of the contour maps of the Quaternary aquifer show that the flow conditions are determined by the Danube main channel as well as the Rivers Iller and Lech. Along the northern border of the Quaternary aquifer in the area near Langenau, a large scale turning of the flow direction in southerly direction indicates larger quantities of infiltrated karst water from the area of the Swabian Alb.

AN US

CR IP T

Fig.3 shows the generated contour maps from the best [a] and worst [d] methods according to the error ranking in Tab.3, as well as from OK [b] and IDW [c] as the most commonly used methods for both aquifers (Jurassic: blue lines; Quaternary: orange lines). In the area where the two aquifers are separated by the low permeable Molasse, the differences of the potential of the two aquifers were calculated to identify possible groundwater exchange zones. Within the blue areas, the Jurassic aquifer has a higher potential (potential rise), within the orange areas, the Quaternary aquifer has a higher potential (potential descent). White areas mean a potential equilibrium.

AC

CE

PT

ED

M

Figure 3: Estimated groundwater contours as well as differences in potential between Jurassic- and Quaternary aquifer based on a: CoKriging D (CoOrdinary Kriging, b: OK (Ordinary Kriging), c: IDW (Inverse Distance Weighting), d: GPI (Global Polynomial Interpolation) and SK (Simple Kriging); U: Ulm; L: Langenau; G: Guenzburg; D: Donauwoerth.

13

ACCEPTED MANUSCRIPT The results show an increase of the karst water potential in the area around Langenau as well as southwest of Donauwoerth (see also Fig. 3). This confirms the results of previous investigations (Villinger 1977; Bierer 1987; Udluft 2000). Further south, the groundwater level of the Quaternary is increasingly above the karst water potential due to the submerging of the Jura layer. The different distances between the observation wells and the inversely proportional weight to the distance on IDW create vast, so-called “bull's eye”, artifacts, which are circular regions of equal values around the known data points. The high gradient thus implied is limited to a small-scale region halfway between the observation wells. In the Jurassic, nearly all methods show an eastward change of direction in the area south of the Danube. Intensity and location of this bending vary greatly within the methods.

PT

ED

M

AN US

CR IP T

Fig. 4 shows a schematic SW-NE cross-section through the unconfined karst of the Swabian Alb. Within the green areas, the maximum distance between the observation wells is less than 2.5 km, yellow less than 5 km, orange less than 10 km. In the cross-section, only the areas with a good observation network, irrespective of the applied method, show plausible results, while in other areas artifacts can be recognized. The best method can not replace an adequate observation network. Here, too, the bullseye effect of IDW is evident. OK and CoOK D show a very similar course in most areas.The course of CoOK D is more disturbed, due to the involvement of the secondary variable, especially the terrain surface. This leads to a probable overfitting especially in areas with high depths to groundwater. The course of SK is approximately equal to the other methods in the southwest area of the cross-section. In the north-east direction, however, the method tends to over- and underestimate the GWL. In the areas of the Kessel valley and the Woernitz valley, all methods fail

CE

Figure 4: Cross-sections depicting topography and groundwater contours generated by Co-Kriging A (Co-Ordinary Kriging), OK (Ordinary Kriging), IDW (Inverse Distance Weighting) and SK (Simple Kriging).

AC

because of the low density of the measuring network which results in an overestimation of the interpolated GWL within the valley. We validated the results additionally with hydrogeological expectations. This includes a comparison between calculated groundwater depth (GWD) and geographic locations of karst springs as well as wetlands and surface waters and finally a comparison of calculated flow accumulation which result from the modeled groundwater surfaces and the locations of real receiving surface streams. Fig. 5 shows the pattern of GWD for the Quaternary aquifer in the lowland fens area Swabian Donaumoos (fen area marked black). Since lowland fens are permanently waterlogged wetlands which are fed by groundwater and rain, the difference between the modeled surface and the DEM surface should ideally be below zero. The yellow patterns show GWD between 3 and 0 m (GWT < DEM), the green patterns show negative GWD between 0 and – 3 m (GWT > DEM). All areas with GWD higher 3m and -3 m were displayed invisible. The soil (Letten) in this area of the Donaumoos were accumulated from clayey lake deposits. These soils form a hydraulic barrier for the groundwater and can lead to confined conditions within the aquifer. Accordingly, negative GWD are not unexpected. CoOK D and OK are able to reproduce the margin of the fen best. For both methods, confined conditions within most areas of the fen were calculated. CoOK D seems to overestimate the GWT at least in the western part of the fen. This is shown

14

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

by the artesian conditions with GWD above 3 m (blank area). These overestimations are also apparent with IDW and to a lesser extend in OK. GPI overestimates the entire southern flank while the north-east flank is underestimated.

15

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 5: Maps of groundwater table depths for the Quaternary in the Donaumoos area ( generated by CoOK A (Co-Ordinary Kriging D), OK (Ordinary Kriging), IDW (Inverse Distance Weighting) and GPI (Global Polynomial Interpolation)

Fig. 6 shows the results for the karst aquifer following the same procedure as for the alluvial aquifer. At a perennial spring,

16

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

the groundwater surface has at least the potential of the terrain. Therefore, the springs can also be used to validate the quality of the interpolated surfaces. Furthermore, the resulting flow direction of the respective methods was calculated with the D8 single-flow algorithm (O'Callaghan & Mark 1984), which is implemented in ArcGIS. For an accurately modeled groundwater surface, the course of the calculated receiving water (shown as purple lines) should be similar to the course of the real receiving waters (shown as light blue lines).

17

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

Figure 6: Maps of groundwater table depths for the Jurassic aquifer and location of perennial karst springs as well as a comparison between rivers and flow directions resulting from the interpolations generated by CoOK D (Co-Ordinary Kriging D), OK (Ordinary Kriging), IDW (Inverse Distance Weighting) and SK (Simple Kriging).

Tab. 4 shows the calculated water exchange of the both aquifers, calculated according to Chapter 2.8. The calculated

18

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

exchange rates within the methods fluctuate of up to a factor of 10. Depending on the methods, an exchange range from 273.96 MCM/yr (Quaternary  Jurassic) to 2.772.63 MCM/yr (Jurassic  Quaternary) could be the consequence. The methods in the upper end of the ranking, however, are in a range between 300 MCM/yr and 500 MCM/yr.

19

ACCEPTED MANUSCRIPT

GPI

RBF

LPI

OK Smooth

OKStandard

SK

CoOK C (Q) /SK (J)

426.67 427.29 -247.12

428.8 425.24 -484.45

427.29 432.69 -974.19

427.6 423.14 -196.97

427.02 427.08 -332.24

426.78 425.72 -227.06

429.02 447.66 -2.451.13

426.57 447.66 -2.772.63

CoOK B

CoOK C

CoOK D

UK

DK

BK

SK (Q)/ LPI (J)

Mean Quaternary GWL [m asl.] Mean Upper Jurassic GWL [m asl.] Exchange btw. Aquifers [MCM/yr]

426.97 428.61 -526.12

426.97 426.87 -429.64

426.46 428.28 -485.76

426.57 427.41 -303.24

427.25 429.61 -640.34

429.38 447.66 -2.462.22

426.78 425.55 -411.84

429.02 423.14 273.96

4. Conclusions

CR IP T

IDW Mean Quaternary GWL [m asl.] Mean Upper Jurassic GWL [m asl.] Exchange btw. Aquifers [MCM/yr]

CoOK A

Table 4: Estimated vertical groundwater exchange rates between the Quaternary alluvial aquifer and the Jurassic karst aquifer. IDW: Inverse Distance Weighting; GPI: Global Polynomial Interpolation; RBF: Radial Basis Functions; LPI: Local Polynomial Interpolation; OK: Ordinary Kriging; SK: Simple Kriging; UK: Universal Kriging; Co-Ordinary Kriging; BK: Empirical Bayesian Kriging. Negative exchange rate: Jurassic  Quaternary

AN US

A total of 9 geospatial and geostatistical (multivariate and univariate) methods were compared with 6 error statistic methods as well as with hydrogeological expectations including the comparison between calculated groundwater depth (GWD) and geographic locations of karst springs as well as wetlands and surface waters. In the multivariate approaches additional surface-elevation data, river levels, and long-term groundwater levels of wells which were not measured at the reference date were included to improve the predictions. The quality of the results for each method has been estimated by qualitative (maps and cross-section) and quantitative (cross-validation) tools.

PT

ED

M

Several important points emerge from the results of the previous section. There is not “a universal superior method” for the interpolations of GWT. Which method performs best is dependent on the number and spatial distribution of the available data points, as well as on the characteristics of the area (type of aquifer, topography etc.). IDW is used most often but is nearly never the “best” method. Especially between remote observation points, IDW tends to create “bull's eye” artifacts in the surface. Geostatistical methods mostly perform better than deterministic methods. In hilly terrain, the estimation of GWT is often problematic because the data set is always sparse in relation to the topographic relief and measuring wells are almost exclusively located in valley region while the groundwater surface is usually a subdued replica of the ground surface elevation (Hoeksema et al. 1989). This often leads to an underestimation of the GWT below mountains ranges and an overestimation of the GWT in valleys. If available, additional data (e.g. topography) can help to improve the results. The consequences of the selection of the interpolation methods can be severe inaccurate or even false computations of recharge, discharge, flow directions, delineation of protections zones etc.). It is advisable to compare the results of different methods, using different error statistics, to compare the results regarding their plausibility, and to give possible ranges for computations based on interpolations results. Of all tested geostatistical methods, it appears that CoOK D (incl. long-term values of monitoring wells not measured at the reference date, DEM data and river levels) is the most accurate method for both investigated aquifers.

AC

CE

Our results show, that the calculated inter-aquifer exchange rates based on different interpolations of groundwater potentials may vary greatly depending on the chosen interpolation method, in our example by a factor of 12.5 or from 273.96 MCM/yr to -2.772.63 MCM/yr comparing the worst case scenarios (highest mean level Quaternary with lowest mean level Jurassic and vice versa). Therefore, the choice of an interpolation method should be made with care, taking into account the different error measures as well as additional data for plausibility control. If a clear choice is not possible, further calculations based on interpolated groundwater contour should at least include an error estimate An inadequate measuring network cannot be replaced by a very good interpolation method. When the data variation is high, the measuring density should be increased to capture the spatial change and therefore to improve the performance of the respective method.

5. Acknowledgements We would like to thank the environmental agencies of Baden-Wuerttemberg and Bavaria (LUBW and LfU) as well the regional water management council of Donauwoerth (WWA-DON) in particularly Dr. Rüdiger Zischak and Cornelius Jakob for providing groundwater level data.

20

ACCEPTED MANUSCRIPT 6. References Aguilar F.J., Agüera F., Aguilar M.A., Carvaja F. (2005): Effects of Terrain Morphology, Sampling Density, and Interpolation Methods on Grid DEM Accuracy. Photogrammetric Engineering & Remote Sensing 71 (7), pp. 805–816. https://dx.doi.org/10.14358/PERS.71.7.805. Ahmadi S.H., Sedghamiz A. (2007): Geostatistical analysis of spatial and temporal variations of groundwater level. Environmental monitoring and assessment 129 (1-3), pp. 277–294. http://dx.doi.org/10.1007/s10661-006-9361-z. Arslan H. (2014): Spatial and temporal distribution of areas with drainage problems as estimated by different interpolation techniques. Water and Environment Journal 28 (2), pp. 203–211. http://dx.doi.org/10.1111/wej.12026.

CR IP T

Bierer S. (1987): 75 Jahre Landeswasserversorgung / Zweckverband Landeswasserversorgung, Stuttgart. 1912 - 1987. Stuttgart. Buhmann M.D. (2003): Radial Basis Functions: Theory and Implementations: Cambridge University Press (Cambridge Monographs on Applied and Computational Mathematics), 272 pp. http://dx.doi.org/10.1017/CBO9780511543241. Chen S., Cowan C.N., Grant P.M. (1991): Orthogonal least squares learning algorithm for radial basis function networks. IEEE transactions on neural networks 2 (2), pp. 302–309. http://dx.doi.org/10.1109/72.80341.

AN US

Chung J.W., Rogers J.D. (2012): Interpolations of groundwater table elevation in dissected uplands. Ground water 50 (4), pp. 598–607. http://dx.doi.org/10.1111/j.1745-6584.2011.00889.x. Cooper H.M., Zhang C., Selch D. (2015): Incorporating uncertainty of groundwater modeling in sea-level rise assessment. A case study in South Florida. Climatic Change 129 (1-2), pp. 281–294. http://dx.doi.org/10.1007/s10584-015-1334-1.

M

Delbari M. (2013): Accounting for exhaustive secondary data into the mapping of water table elevation. Arabian Journal of Geoscience 7 (10), pp. 4221–4233. http://dx.doi.org/10.1007/s12517-013-0986-2. Flinspach D., Haakh F., Locher A., (1997): Das württembergische Donauried - seine Bedeutung für die Wasserversorgung, Landwirtschaft und Naturschutz, 271 p. Zweckverband Landeswasserversorgung, Stuttgart.

ED

Guekie A.T., Marache A., Lastennet R., Breysse D. (2016): Geostatistical investigations for suitable mapping of the water table. The Bordeaux case (France). Hydrogeology Journal 24 (1), pp. 231–248. http://dx.doi.org/10.1007/s10040-015-1316-4.

PT

Hengl T. (2009): A practical guide to geostatistical mapping. JRC Scientific and Technical Reports, 271 p. 2nd extended edition. Amsterdam.

CE

Hoeksema R.J., Clapp R.B., Thomas A.L., Hunley A.E., Farrow N.D., Dearstone K.C. (1989): CoKriging model for estimation of water table elevation. Water Resources Research. 25 (3), pp. 429–438. http://dx.doi.org/10.1029/WR025i003p00429. Hua Z., Dehai M., Cheng W. (2009): Optimization of the Spatial Interpolation for Groundwater Depth in Shule River Basin, DBLP Conference: 2009 Intern. Conference on Environmental Science Application Technology, pp. 415–418. http://dx.doi.org/10.1109/ESIAT.2009.478.

AC

Johnston K. (2004): ArcGis 9. Using ArcGis geostatistical analyst. Redlands, California. Esri Press (GIS by ESRI). Li J., Heap A.D. (2008): A review of spatial interpolation methods for environmental scientists. Canberra. Geoscience Australia Record 2008/23, 137 pp. Li J., Heap A.D. (2011): A review of comparative studies of spatial interpolation methods in environmental sciences. Performance and impact factors. Ecological Informatics 6 (3-4), pp. 228–241. http://dx.doi.org/10.1016/j.ecoinf.2010.12.003. Möhler F., Dinse S., Hermsdorf A. (2014): Grundwassergleichenplan für Brandenburg – Interpolation mittels Kriging mit externer Drift. Grundwasser 19 (3), pp.189–199. http://dx.doi.org/10.1007/s00767-014-0255-7. O'Callaghan J. F., Mark D.M. (1984): The extraction of drainage networks from digital elevation data. Computer Vision. Graphics, & Image Processing 28 (3), pp. 323–344. http://dx.doi.org/10.1016/S0734-189X(84)80011-0

21

ACCEPTED MANUSCRIPT Oliver M.A., Webster R., McGrath S.P. (1996): Disjunctive Kriging for Environmental Management. Environmetrics 7 (3), pp. 333–357. http://dx.doi.org/10.1002/(SICI)1099-095X(199605)7:3<333::AID-ENV209>3.0.CO;2-V. Rivoirard J. (1994): Introduction to disjunctive Kriging and non-linear geostatistics. Oxford: Oxford Univ. Press (Spatial Information Systems), 180 pp. Sadat Noori S.M., Ebrahimi K., Liaghat A.M., Hoorfar A.H. (2012): Comparison of different geostatistical methods to estimate groundwater level at different climatic periods. Water Environmental Journal 27 (1), pp. 10–19. http://dx.doi.org/10.1111/j.1747-6593.2012.00321.x. Snyder D.T. (2008): Estimated Depth to Ground Water and Configuration of the Water Table in the Portland, Oregon Area. Scientific Investigations Report 2008–5059. 40 pp. http://pubs.usgs.gov/sir/2008/5059/.

CR IP T

Sun Y., Kang S., Li F., Zhang L (2009): Comparison of interpolation methods for depth to groundwater and its temporal and spatial variations in the Minqin oasis of northwest China. Environmental Modelling & Software 24 (10), pp 1163–1170. http://dx.doi.org/10.1016/j.envsoft.2009.03.009. Tapoglou E., Karatzas G.P., Trichakis I.C., Varouchakis E.A. (2014): A spatio-temporal hybrid neural network-Kriging model for groundwater level simulation. Journal of Hydrology 519, pp. 3193–3203. http://dx.doi.org/10.1016/j.jhydrol.2014.10.040.

AN US

Udluft P. (2000): Das Grundwasser im schwäbischen Donauteil. Hydrologisch-hydrogeologische Untersuchung mit Erstellung eines Grundwassermodells im Maßstab 1:25000/50000 im Donautal zwischen Ulm/Neu-Ulm und Neuburg an der Donau. München: Bayer. Industrieverb. Steine und Erden Fachabt. Sand- und Kiessteine (Schriftenreihe der bayerischen Sand- und Kiesindustrie, H. 11), 102 pp. USGS (2014): Shuttle Radar Topography Mission (SRTM) 1 Arc-Second. N48E009V3 - N48E010V3. U.S. Geological Survey data. https://earthexplorer.usgs.gov/. Varouchakis E.A., Hristopulos D.T. (2013): Comparison of stochastic and deterministic methods for mapping groundwater level spatial variability in sparsely monitored basins. Environmental monitoring and assessment 185 (1), pp. 1–19. http://dx.doi.org/10.1007/s10661-012-2527-y.

M

Villinger E. (1977): Über Potentialverteilung und Strömungssysteme im Karstwasser der Schwäbischen Alb. (Oberer Jura, SW-Deutschland). Stuttgart: Schweizerbart (Geologisches Jahrbuch : Reihe C, Hydrogeologie, Ingenieurgeologie, H. 18), 96 pp.

ED

Wang S., Huang G.H., Lin Q.G., Li Z. Zhang H., Fan Y.R. (2014): Comparison of interpolation methods for estimating spatial distribution of precipitation in Ontario, Canada. International Journal of Climatology 34 (14), pp. 3745–3751. http://dx.doi.org/10.1002/joc.3941.

PT

Winkler H.A.G. (1972): Das Grundwasser im Nördlinger Ries unter Berücksichtigung der hydrologischen und hydrochemischen Beziehungen zum Speichergestein. Dissertation. Technische Universität München, München. Fakultät f. Allg. Wissenschaften, 202 pp.

CE

Xiao Y., Gu X., Yin S., Shao J., Cui Y., Zhang Q., Niu Y. (2016): Geostatistical interpolation model selection based on ArcGIS and spatio-temporal variability analysis of groundwater level in piedmont plains, northwest China. SpringerPlus 5, pp. 1-15. http://dx.doi.org/10.1186/s40064-016-2073-0.

AC

Xie Y., Chen T.B., Lei M., Yang J., Guo Q.J., Song B., Zhou X.Y. (2011): Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: accuracy and uncertainty analysis. Chemosphere 82 (3), pp. 468-476. http://dx.doi.org/10.1016/j.chemosphere.2010.09.053. Yao L., Huo Z., Feng S., Mao X., Kang S., Chen J., (2013): Evaluation of spatial interpolation methods for groundwater level in an arid inland oasis, northwest China. Environmental Earth Science 71 (4), pp. 1911-1924. http://dx.doi.org/10.1007/s12665-013-2595-5. Zimmerman D., Pavlik C., Ruggles A., Armstrong M.P. (1999): An Experimental Comparison of Ordinary and Universal Kriging and Inverse Distance Weighting. In: Mathematical Geology 31 (4), pp. 375-390. http://dx.doi.org/10.1023/A:1007586507433.

22

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Graphical abstract

23