Geostatistical analysis and space-time models of aquifer levels: application to Mires hydrological basin in the prefecture of Crete

Abstract

The accurate representation of groundwater levels in an aquifer is very important for groundwater modeling and effective groundwater resources management. However, the number and spatial distribution of monitoring sites in a given aquifer are not always sufficient to accurately represent the water table. Predictions of groundwater level at unvisited locations of an aquifer can be obtained by applying geostatistical methods on the available groundwater level data and thus the free surface of the aquifer can be reliably mapped. In sparsely monitored basins, accurate mapping of the spatial variability of groundwater level requires the interpolation of scattered data. This thesis aims to present new modeling tools that help to better monitor and predict the groundwater level in sparsely gauged basins. The specific area of focus is the Mires basin of the Mesara valley in the island of Crete (Greece). The study area is a sparsely sampled basin that has limited groundwater resources which are vital for the area’s welfare; spatiotemporal variations of groundwater level are important for developing management and monitoring strategies. Efficient groundwater management in the basin is crucial in light of regional climate change model estimates showing a substantial risk of desertification for Crete. Our goal is to construct accurate spatial and spatiotemporal models of the basin’s groundwater level. Therefore, spatial and spatiotemporal models for the accurate representation of the groundwater level variability in already vulnerable areas with low groundwater resources, like Mires basin, need to be developed in order to identify the susceptible locations, to estimate the groundwater level distribution spatially and spatiotemporally and to provide input for potential groundwater resources management plans. The main data used in this research consist of seventy hydraulic head measurements (wet period of 2002-2003 hydrological year) which are unevenly distributed and mostly concentrated along a temporary river and time series performance consisting of biannual groundwater level data from tenboreholes (1981-2003). After the year 2003 observations of a shorter number of wells are available biannually. This thesis initially presents a comparison of deterministic interpolation methods, i.e., Inverse Distance Weight (IDW) and Minimum Curvature (MC), with stochastic methods, i.e., Ordinary Kriging (OK), Universal Kriging (UK) and Kriging with Delaunay Triangulation (DK). We evaluate the performance of the interpolation methods with respect to different statistical cross validation measures. The Spartan variogram family is implemented for the first time to hydrological data and is shown to be optimal with respect to the stochastic interpolation methods such as OK, UK, DK and Residual Kriging (RK) applied in this dataset. The three stochastic methods (OK, DK, UK) perform overall better than the deterministic counterparts (IDW, MC). DK, which is herein for the first time applied to hydrological data, yields the most accurate cross validation estimate for the lowest value in the dataset. OK and UK lead to smooth isolevel contours, while DK and IDW generate more edges. The stochastic methods also deliver estimates of prediction uncertainty. The present research study also investigates the application of non-linear normalizing data transformations in conjunction with Ordinary Kriging (OK) for the accurate prediction of groundwater level spatial variability in a sparsely gauged basin. We investigate three established normalizing methods, Gaussian Anamorphosis, Trans-Gaussian Kriging and Box-Cox to improve the estimation accuracy. The first two are for the first time applied to groundwater level data. All three methods improve the mean absolute prediction error compared with the application of OK to the non-transformed data. In addition, a Modified Box-Cox (MBC) transformation is proposed and applied to normalize the hydraulic heads. MBC in conjunction with OK is found to be the optimal spatial model based on leave-one-out cross-validation. The recently established Spartan semivariogram family provides the optimal model fit to the transformed data. Trend functions, as previous studies have shown, improve the accuracy of interpolation. Therefore, we propose that the prediction of the hydraulic head spatial variability in Mires basin can be improved by incorporating in the trend function local properties. Firstly the distance of the prediction points from the temporary river crossing the basin is incorporated in the trend function and secondly a component based on the generalized Thiem’s equation for multiple wells. Residual Kriging is performed based on these two spatial trend models as well as using the novel MBC transformation to normalize the residuals and the flexible Spartan semivariogram family to optimally determine their spatial correlation. Both proposed spatial models’ improve significantly the cross validation measures compared to the other Krigingbased methods tested. We also present maps of the groundwater level spatial variability and the estimation variance in Mires basin obtained by means of the optimal spatial models. We determine that the optimal spatial interpolation approach is based on Residual Kriging using the non-differentiable Spartan semivariogram model for the normalized (MBC) fluctuations. Our proposal is supported by the results of cross validation analysis. The suggested methodology is applicable to other unconfined aquifers as well. The non-differentiable property of the Spartan semivariogram model is interpreted herein as the result of a deposition-removal process that leads to a fractional Brownian motion, fBm-like, behavior of the groundwater level surface. The overall goal of this dissertation is to use stochastic methods for the monitoring and prediction of the groundwater level in space and time. Following the spatial interpolation, first we model the temporal variation of groundwater level with a discrete time autoregressive exogenous variable model (ARX) model. In this study pumping data are used in addition to precipitation measurements. The ARX model is embedded in a discrete-time Kalman filter to estimate the model parameters andpredict the optimal mean annual groundwater level. The ARX model is calibrated for the years 1981 to 2006 and is then used to predict the mean annual groundwater level in the basin for recent years (2007-2010). The predictions are validated with the available annual averages reported by the local authorities. Secondly, we use a spatiotemporal geostatistical analysis of the groundwater level using space-time Residual Kriging (STRK). The space-time trend is calculated using the product function of the estimated temporal trend from a weighted moving average filter and the spatial trend determined from the closest distance of the measurement locations from the river bed. A space-time experimental semivariogram is determined from the biannual (wet and dry period) groundwater level fluctuations time series between the years 1981 and 2003 at the ten sampling stations. We model the semivariogram with separable and non-separable theoretical spatiotemporalsemivariogram functions. STRK is used to predict the groundwater level for selected hydrological periods at each sampling station in the time period (2004-2010) biannually. Maps of groundwater level predictions and of prediction accuracy are desirable and significant in order to assess the groundwater level spatiotemporal variability, whether observed changes in water-table levels are statistically significant and finally to identify additional locations where further monitoring is needed to increase the accuracy of the maps. All the methodologies and tools presented in this thesis are implemented by original code developed by the author and run in Matlab® programming environment (Matlab v.7.5 on Microsoft Windows XP).

Η ακριβής εκτίμηση της στάθμης του υδροφόρου ορίζοντα αποτελεί σημαντικό παράγοντα για ένα ολοκληρωμένο σχέδιο διαχείρισης των υδατικών πόρων μιας υδρολογικής λεκάνης. Η στάθμη των υπόγειων υδάτων σε έναν υδροφορέα ελέγχεται συνήθως με βάση το υδραυλικό ύψος στις θέσειςεγκατεστημένων πιεζόμετρων και ενεργών γεωτρήσεων σε κατάσταση ηρεμίας. Σε πολλές περιπτώσεις δεν παρακολουθούνται όλες οι γεωτρήσεις λόγω κόστους ή παράλειψης των αρμόδιων φορέων, ενώ σε άλλες περιπτώσεις ο αριθμός των γεωτρήσεων και των πιεζομέτρων είναι ανεπαρκής. Επομένως ο αριθμός των διαθέσιμων μετρήσεων είναι συχνά μικρός και αραιά κατανεμημένος στην περιοχή μελέτης με αποτέλεσμα να μην αντιπροσωπεύεται επαρκώς η στάθμη του υδροφορέα. Εκτιμήσεις της στάθμης σε θέσεις όπου δεν υπάρχουν παρατηρήσεις μπορούν ναπραγματοποιηθούν με την εφαρμογή γεωστατιστικών μεθόδων στα διαθέσιμα δεδομένα, προκειμένου να χαρτογραφηθεί με ακρίβεια ο υδροφόρος ορίζοντας του υδροφορέα. Βοηθητικές πληροφορίες χρησιμοποιούνται συχνά για να ενισχύσουν τις εκτιμήσεις της στάθμης των υπόγειων νερών, όπως η τάση της φυσικής μεταβολής του επιπέδου των υπόγειων νερών (η οποία προσεγγίζεται συνήθως από πολυώνυμα με βάση τις χωρικές συντεταγμένες), η βροχόπτωση και το υψόμετρο. Η χρήση βοηθητικών χωρικών μεταβλητών έχει αποδειχθεί ότι βελτιώνει τις εκτιμήσειςυδραυλικών υψών. Το νησί της Κρήτης διαθέτει οριακούς υπόγειους υδατικούς πόρους, οι οποίοιχρησιμοποιούνται εκτενώς για γεωργικές δραστηριότητες και ύδρευση. Η κοιλάδα της Μεσαράς, η οποία βρίσκεται στο νότιο τμήμα του περιφερειακού διαμερίσματος Ηρακλείου και καλύπτει μια έκταση 398 km2 , είναι η μεγαλύτερη και παραγωγικότερη κοιλάδα του νησιού. Η υπερεκμετάλλευση κατά τη διάρκεια των προηγούμενων τριάντα ετών έχει οδηγήσει σε μια δραματική μείωση, πάνω από τριάντα πέντε μέτρα, στη στάθμη των υπόγειων νερών. Οι πιθανές μελλοντικές κλιματολογικές αλλαγές στην περιοχή της Μεσογείου, τα σενάρια πιθανής ερημοποίησης και η εκτενής γεωργική δραστηριότητα προκαλούν έντονη ανησυχία σχετικά με την αειφορία των υδατικών πόρων της περιοχής. Η παρούσα διατριβή εστιάζει στην υδρολογική λεκάνη Μοιρών της κοιλάδας της Μεσσαράς για δύο λόγους: α) διαθεσιμότητα υδρογεωλογικών και υδρολογικών δεδομένων και β) διότι αποτελείται κυρίως από προσχωματικές αποθέσεις κάτι που εξασφαλίζει σε μεγάλο βαθμό υδρογεωλογική ομοιογένεια. Η παρούσα διατριβή έχει στόχο να παρουσιάσει καινοτόμες μεθοδολογίες χωρικής και χρονικής γεωστατιστικής ανάλυσης αλλά και μεθόδους–εργαλεία αποσκοπώντας στην εκτίμηση και χαρτογράφηση της χωρικής και χρονικής μεταβλητότητας της στάθμης των υπόγειων νερών της λεκάνης με τη βέλτιστη δυνατή ακρίβεια λαμβάνοντας υπόψιν και την αβεβαιότητα των εκτιμήσεων. Ο κύριος στόχος της έρευνας είναι να αναπτυχθούν χωρικά μοντέλα για τον ακριβή προσδιορισμό της στάθμης, ο υπολογισμός των αντίστοιχων ισοδυναμικών καμπυλών και ο προσδιορισμός των πιο ευάλωτων περιοχών του υδροφορέα. Σε δεύτερο στάδιο, τα μοντέλα επεκτείνονται σε δυναμικές (space-time) καταστάσεις για τη μοντελοποίηση της χωροχρονικής μεταβολής και τη δυνατότητα εκτίμηση της μελλοντικής στάθμης. Η χωρική εξάρτηση τόσο των δεδομένων στάθμης όσο και των διακυμάνσεων που προκύπτουν από τα μοντέλα τάσης που αναπτύσσονται στην παρούσα διδακτορική διατριβή μελετώνται με την προσαρμογή εμπειρικών ημιβαριογραμμάτων σε γνωστά κλασικά πρότυπα συναρτήσεων, στο πρότυπο συναρτήσεων Matérn και στο μοντέλο Σπαρτιάτικου ημιβαριογράμματος το οποίο χρησιμοποιείται για πρώτη φορά σε υδρολογικά δεδομένα. Τα μοντέλα Σπαρτιάτικουτύπου και Matérn περιλαμβάνουν περισσότερες παραμέτρους σε σχέση με τα κλασικά πρότυπα, γεγονός που ευνοεί στη βέλτιστη προσαρμογή στο πειραματικό ημιβαριόγραμμα. Το Σπαρτιάτικο ημιβαριόγραμμα και η αντίστοιχη συνάρτηση συνδιασποράς έχουν αναπτυχθεί πρόσφατα και αποτελούν μια επιτυχώς εναλλακτική πρόταση στον προσδιορισμό της χωρικής εξάρτησης των δεδομένων. Αρχικά αξιολογείται η χρήση γνωστών και ευρέως εφαρμοσμένων μεθόδων χωρικής παρεμβολής όπως το κανονικό Kriging (ordinary Kriging-ΟΚ), το γενικευμένο Kriging (universal Kriging-UK), η μέθοδος σταθμισμένων αντίστροφων αποστάσεων (inverse distance weight-IDW), η μέθοδος ελάχιστης καμπυλότητας (minimum curvature-MC) και η μέθοδος Kriging με τριγωνοποίηση κατά Delaunay (Kriging with Delaunay Triangulation-DK). Η τελευταία διερευνάται για πρώτη φοράσε υδρολογικά δεδομένα. Μάλιστα εξετάζεται και παρουσιάζεται η εφαρμογή της μεθόδου χρησιμοποιώντας δύο διαφορετικές γειτονίες εκτίμησης που προκύπτουν από τον τρόπο επιλογής των γειτονικών τιμών του εκτιμώμενου σημείου με βάση την πρωτεύουσα και δευτερεύουσα τριγωνοποίηση των γειτονικών τιμών. Η σύγκριση των συγκεκριμένων πέντε μεθόδων στο ίδιο δείγμα πραγματοποιείται για πρώτη φορά σύμφωνα με τη μέχρι πρόσφατη βιβλιογραφία. Τα αποτελέσματα της σύγκρισης υποδεικνύουν ότι οι στοχαστικές μέθοδοι (OK, DK, UK) λειτουργούν με μεγαλύτερη αποτελεσματικότητα σε σχέση με τις αντίστοιχες προσδιοριστικές (IDW, MC). Επίσης το Σπαρτιάτικο ημιβαριόγραμμα αποδεικνύεται με βάση διαφορετικά στατιστικά μέτρα το βέλτιστο για την εφαρμογή της γωστατιστικής εκτίμησης. Οι μέθοδοι OK και UK οδηγούν σε πιο ομαλές ισοδυναμικές σε αντίθεση με τις μεθόδους DK και IDW. Η μέθοδος DK η οποία χρησιμοποιείται για πρώτη φορά σε υδρολογικά δεδομένα υπολογίζει ακριβέστερα την ελάχιστη τιμή του δείγματος βάση της μεθόδου διασταυρωμένης επιβεβαίωσης. Οι στοχαστικές μέθοδοι πλεονεκτούν σεσχέση με τις προσδιοριστικές καθώς υπολογίζουν την αβεβαιότητα των εκτιμήσεων η οποία και μπορεί να υποδηλώσει περιοχές της υπό μελέτη λεκάνης όπου περισσότερες μετρήσεις απαιτούνται για γεωστατιστική ανάλυση μεγαλύτερης ακρίβειας. Η προκαταρκτική γεωστατιστική ανάλυση των δεδομένων υδραυλικών υψών έδειξε ότι αυτά δεν ακολουθούν την κανονική κατανομή ωστόσο δεν είναι απαγορευτική η χρήση τους σε γραμμική γεωστατιστική ανάλυση. Για να δημιουργηθεί όμως ένα χωρικό μοντέλο με τη βέλτιστη δυνατή αποτελεσματικότητα για τα υδραυλικά ύψη εξετάζονται διάφορες μη γραμμικές προσεγγίσεις κανονικοποίησης των δεδομένων σε συνδυασμό με τη μέθοδο εκτίμησης kriging. Η μέθοδος Box-Cox και η μέθοδος της Γκαουσιανής Αναμόρφωσης (GaussianAnamorphosis-GA) χρησιμοποιούνται για το μετασχηματισμό των δεδομένων στην κανονική (γκαουσιανή) κατανομή πιθανότητας. Προτείνεται επίσης μια νέα τροποποιημένη εκδοχή της μεθόδου Box-Cox (modified Box-Cox) η οποία βασίζεται στο συντελεστή κύρτωσης και ασυμμετρίας της παρατηρούμενης κατανομής. Χρησιμοποιείται ακόμη η μέθοδος Trans-Gaussian Kriging η οποία ενσωματώνει τη συνάρτηση μετασχηματισμού των δεδομένων στον εκτιμητή της χωρικήςπαρεμβολής. Η χρήση των εν λόγω μη γραμμικών μοντέλων αποτελεί καινοτομία στην υδρολογία. Η μέθοδος modified Box-Cox σε συνδυασμό με τη μέθοδο kriging και το Σπαρτιάτικο ημιβαριόγραμμα αποτελούν το βέλτιστο χωρικό μοντέλο εκτίμησης της στάθμης του υδροφορέα βάση σύγκρισης με τις υπόλοιπες μεθόδους. Παρουσιάζονται χάρτες ισοδυναμικών καμπυλών και αβεβαιότητας εκτιμήσεων με όλες τις στοχαστικές μεθόδους που εξετάστηκαν. Επίσης προτείνεται και διερευνάται ένα καινοτόμο χωρικό μοντέλο τάσης για τη στάθμη των υπογείων υδάτων που βασίζεται στα τοπικά χαρακτηριστικά της λεκάνης των Μοιρών, και συγκεκριμένα στην απόσταση των γεωτρήσεων από την κοίτη του ποταμού. Το μοντέλο επίσης ενσωματώνει το αντίστοιχο υψόμετρο στη θέση της κάθε γεώτρησης. Μια δεύτερη καινοτόμα πρόταση χωρικού μοντέλου τάσης αποτελεί η χρήση των πολλαπλών πηγαδιών άντλησης στον υδροφορέα των Μοιρών. Η εφαρμογή της εξίσωσης Thiem για πολλαπλά πηγάδια άντλησης μπορεί να αποδώσει την τάση της στάθμης της λεκάνης. Η χρήση των προτεινόμενων μεθόδων για τον υπολογισμό της τάσης της στάθμης του υδροφορέα σε συνδυασμό με τη μέθοδο του υπολειμματικού kriging (Residual Kriging) οδηγεί στη βέλτιστη εκτίμηση του υδραυλικού ύψους σε σημεία της λεκάνης όπου δεν υπάρχουν παρατηρήσεις. Η μέθοδος συνοδεύεται από τη χρήση της καινοτόμου μεθόδου κανονικοποίηση δεδομένων modified Box-Cox (για τις διακυμάνσεις) και του μη διαφορίσιμου Σπαρτιάτικου ημιβαριόγραμματος για τον προσδιορισμό της χωρικής εξάρτησης των διακυμάνσεων. Οι προτάσεις μας υποστηρίζονται από τα αποτελέσματα διασταυρωμένης επιβεβαίωσης. Το Σπαρτιάτικο ημιβαριόγραμμα αποτελεί το βέλτιστο μοντέλο προσδιορισμού της χωρικής εξάρτησης των δεδομένων για κάθε μεθοδολογία που εξετάστηκε στην παρούσα διατριβή. Ωστόσο η συνάρτηση αυτή είναι μη διαφορίσιμη. Η ιδιότητα αυτή ερμηνεύεται ως το αποτέλεσμα μιας διαδικασίας εμπλουτισμού και άντλησης του υδροφορέα η οποία οδηγεί τον υδροφόρο ορίζοντα σε συμπεριφορά κλασματικής κίνησης Brown (fractional Brownian motion). Η μοντελοποίηση της χρονικής μεταβολής της στάθμης του υδροφορέαπραγματοποιείται με τη χρήση ενός μοντέλου αυτοσυσχέτισης το οποίο ενσωματώνει εξωγενή πληροφορία από μεταβλητές όπως η βροχόπτωση, η παροχή αντλήσεων και η εξατμισοδιαπνοή. Το συγκεκριμένο μοντέλο έχει χρησιμοποιηθεί στην αρχική του μορφή χωρίς την παροχή αντλήσεων. Στην παρούσα διατριβή διερευνάται η χρήση και της μεταβλητής αυτής για πρώτη φορά με επιτυχή αποτελέσματα. Το μοντέλο αυτοπαλινδρόμησης ενσωματώνεται σε ένα διακριτό χρονικά φίλτρο Kalman για την εκτίμηση των παραμέτρων αλλά και για την πρόβλεψη της βέλτιστης στάθμης τουυδροφορέα. Τα κύρια στοιχεία που χρησιμοποιούνται στην παρούσα διατριβή αποτελούνται από εβδομήντα (70) μετρήσεις στάθμης που αφορούν την υγρή περίοδο του υδρολογικού έτους 2002-2003 και από στάθμες 10 γεωτρήσεων που καλύπτουν το χρονικό διάστημα 1981-2003 σε υγρή και ξηρή περίοδο. Οι μετρήσεις κατανέμονται ανομοιόμορφα και εστιάζονται κατά μήκος του ποταμού Γεροπόταμου που διασχίζει τη λεκάνη. Από το 2003 και μετά, μετρήσεις πραγματοποιήθηκαν μόνοσε μικρό αριθμό επιλεγμένων γεωτρήσεων ενώ από το 2003 δύο τηλεμετρικοί σταθμοί λειτουργούν στην περιοχή παρακολουθώντας τη μεταβολή της στάθμης των υπογείων υδάτων. Η ακριβής χωρική εκτίμηση της στάθμης σε συνδυασμό με τη μοντελοποίηση της χρονικής μεταβολής δημιουργούν τις συνθήκες για ένα ολοκληρωμένο χωροχρονικό μοντέλο το οποίο θα περιγράφει με ακρίβεια τη στάθμη του υδροφορέα και θα μπορεί να εκτιμά και τη μελλοντική συμπεριφορά του. Για το σκοπό αυτό οι τάσεις της χρονικής μεταβολής της στάθμης μελετώνται και προσδιορίζονται από έναμοντέλο σταθμισμένου κινούμενου μέσου όρου (weighted moving average). Η χωροχρονική τάση προσδιορίζεται από το συνδυασμό του σταθμισμένου κινούμενου μέσου όρου και της απόσταση των γεωτρήσεων από την κοίτη του ποταμού. Στη συνέχεια η χωροχρονικής εξάρτηση των διακυμάνσεων των μετρήσεων υπολογίζεται με τη βοήθεια του χωροχρονικού εμπειρικού ημιβαριογράμματος. Η μοντελοποίηση του πραγματοποιείται με τη χρήση διαχωριζόμενων και μη διαχωριζόμενων χωροχρονικών συναρτήσεων. Η επαλήθευση εκτίμηση και πρόβλεψη στάθμηςπραγματοποιείται με την εφαρμογή του χωροχρονικού υπολειμματικού kriging (Residual Kriging). Η χρήση μη διαχωρίσιμων χωροχρονικών συναρτήσεων αποτελεί επίσης καινοτομία στην υδρολογία.Οι μεθοδολογίες που αναφέρονται παραπάνω δύναται να εφαρμοστούν και σε άλλες υδρολογικές λεκάνες με παρόμοια χαρακτηριστικά όπως αυτά της λεκάνης των Μοιρών. Η υλοποίηση των προτεινόμενων μεθόδων πραγματοποιήθηκε από το συγγραφέα σε πρωτότυπο κώδικα στο προγραμματιστικό περιβάλλον Matlab® ενώ η αξιολόγηση των γεωστατιστικών μεθόδων πραγματοποιείται με τη μέθοδο της διασταυρωμένης επιβεβαίωσης χρησιμοποιώντας διάφορα στατιστικά μέτρα επίδοσης. Πιστεύουμε ότι η παρούσα διδακτορική διατριβή συνεισφέρει στηνεπιστήμη της γεωστατιστικής αλλά και στην υδρολογική διερεύνηση της περιοχής μελέτης. Τα αποτελέσματα και οι μέθοδοι που παρουσιάζονται μπορούν να χρησιμοποιηθούν σε συνδυασμό με άλλες γεωεπιστήμες για την ολοκληρωμένη διαχείριση των υδατικών πόρων της λεκάνης των Μοιρών και της ευρύτερης περιοχής της Μεσσαράς.