Ανάλυση και διερεύνηση της βιώσιμης διαχείρισης των υπόγειων νερών με μαθηματικά μοντέλα σε περιβάλλον GIS

Abstract

The objective of this dissertation is a critical review of all the required stages for the development of a sound groundwater mathematical model. Towards this objective a groundwater model is developed and implemented for the management of the Pieria aquifer, Greece. The stages of the modeling process which are analyzed in this dissertation are the following: - development of the conceptual model, - construction of the numerical model, which includes the definition of the space and time discretization, the initial and boundary conditions, the recharge and discharge of the system, the representation of hydrogeologic features of the system, the spatial structure of the aquifer’s hydrogeologic parameters, the selection of the method for the solution of the groundwater flow equation and the convergence criteria for this method, - calibration of the model, - evaluation of the model, - sensitivity analysis, - validation of the calibrated model and – application of the model to predict the system’s response under various conditions. The model which was developed for the Pieria aquifer is a two-dimensional transient flow model. It covers an area of 256 km² and it has 1759 active cells, with a size of 350 m x 350 m each. Software packages which were used for the construction of the model are: MODFLOW 2000, Argus ONE, ArcGIS and PEST. The zonation method was used for the description of the hydrogeologic parameters’ spatial variability. The grid was divided in 11 zones of uniform hydraulic conductivity in the x-direction, 6 zones of uniform horizontal anisotropy, 4 zones of uniform specific storage and 4 zones of uniform specific yield. The values of the 25 hydrogeologic parameters resulted from the calibration process. The model was calibrated against 1173 piezometric head observations and 26 spring flow observations with the Gauss-Levenberg-Marquardt method. The piezometric head and the spring flow observations were collected from September 1992 to April 1995, on a monthly basis. Both qualitative and quantitative criteria were used for the evaluation of the calibrated model. The results which we obtained show that the model represents adequately the simulated system. The effect of each one of the 38 model parameters (25 hydrogeologic and 13 recharge parameters) was examined with the sensitivity analysis process. A second set of monthly data (118 head observations and 4 spring discharge observations, from May 1995 to August 1995 were used for the validation of the calibrated model. The validation process verified that the calibrated model describes accurately the water dynamics of the Pieria aquifer. Finally, the calibrated model was applied for a period of 3 years to assess the consequences of two pumping scenarios concerning the location of irrigation wells. The proper time discretization was examined by using a number of uniform and geometrically increased time steps. The comparison of total execution time and computed heads revealed that increasing the length of the time step geometrically, can reduce execution time significantly, but the corresponding heads diverge from those computed with uniform time steps. It is concluded that the time discretization depends on the special characteristics of the model and must be investigated thoroughly in any groundwater mathematical model. The initial condition of the computations must be also examined very carefully. An arbitrarily defined initial condition (i.e. interpolation from field measurements), sometimes can influence the results until the end of the whole simulation period. The known head boundary condition constitutes an inexhaustible water source for the model, so its contribution to the model’s water balance should also be examined and taken into account during the calibration of the model. Sensitivity analysis of this boundary condition, in order to assess the magnitude and the spatial distribution of its influence, is an essential prerequisite for the model construction. The particularity of this boundary condition should also be taken into consideration during predictive simulations. Plots of the computed recharge/discharge from every known head boundary cell with respect to time can be a valuable tool for the estimation of the known head cell values in predictive simulations. It is concluded that the integration of this boundary condition in the calibration process leads to a more accurate representation of the real system. In our case this was accomplished using artificial observations of zero inflow from the known head boundary cells. The head convergence criterion is very important and depends on the solution method. It must always be examined in conjunction with the computed water balance error. In our case this error never exceeded 0.1%. It was found that the use of different kind of observations in the calibration process decreased the correlation between the estimated parameters and thus it led to a unique solution of the inverse problem. Model efficacy in predicting the response of the system depends on the quality and quantity of available data. A transient model should be calibrated against a considerable number of observations with satisfactory spatial and temporal distribution, representing the system response under various conditions. The evaluation of a calibrated model must be based on statistical criteria referring to the residuals between observations and corresponding computed values, to the absolute residuals and to the final parameter values resulted from the calibration process. The water balance of the model should also be considered. Most of the commonly used statistical criteria come from time series analyses where large data sets are usually available. In groundwater models, however, availability of such data sets is rare. It is thus important to use several criteria to get a better evaluation of the model. Software’s automated procedures for data handling and parameter default values often mislead the development of a groundwater model. Appropriate interpolation techniques and parameter values are problem depended and should be decided after investigation.

Αντικείμενο της διατριβής αυτής είναι η ανάλυση και διερεύνηση όλων των σταδίων της διαδικασίας κατάρτισης και εφαρμογής ενός μαθηματικού μοντέλου υπόγειων νερών. Μέρος αυτής της ανάλυσης αποτελεί και ο εντοπισμός και η αντιμετώπιση σφαλμάτων που συχνά πραγματοποιούνται και οδηγούν σε τελείως λανθασμένα αποτελέσματα. Η διερεύνηση πραγματοποιείται μέσω της κατάρτισης κι εφαρμογής ενός μοντέλου για τα υπόγεια νερά της πεδιάδας Πιερίας. Στα πλαίσια της διατριβής αναλύθηκαν τα στάδια: - ανάπτυξης του θεμελιώδους μοντέλου, - κατάρτισης του μαθηματικού μοντέλου, το οποίο περιλαμβάνει: την επιλογή της χωρικής και χρονικής διακριτοποίησης, τον ορισμό της αρχικής συνθήκης, των οριακών συνθηκών, του εμπλουτισμού του συστήματος, της εκφόρτισής του, του τρόπου προσομοίωσης επιμέρους υδρογεωλογικών χαρακτηριστικών του συστήματος, την περιγραφή της χωρικής δομής των υδρογεωλογικών παραμέτρων, την επιλογή της μεθόδου επίλυσης της εξίσωσης που περιγράφει την υπόγεια ροή και των κατάλληλων κριτηρίων σύγκλισης της μεθόδου, - ρύθμισης του μοντέλου, - ανάλυσης ευαισθησίας των παραμέτρων του, - επαλήθευσής του και τέλος – εφαρμογής του για την πρόβλεψη της αντίδρασης του συστήματος σε διάφορες συνθήκες λειτουργίας του. Το μοντέλο που αναπτύχθηκε καλύπτει μια έκταση περίπου 256 km², προσομοιώνει διδιάστατη μη μόνιμη ροή και αποτελείται από 1759 ενεργά κελιά, διαστάσεων 350m x 350m. Για την κατάρτισή του χρησιμοποιήθηκε το MODFLOW 2000, μέσα από το γραφικό περιβάλλον του Argus-ONE, σε συνεργασία με το σύστημα γεωγραφικής πληροφορίας του ArcGIS. Για την περιγραφή της χωρικής δομής των υδρογεωλογικών παραμέτρων του υδροφορέα, χρησιμοποιήθηκε η μέθοδος διαχωρισμού κάθε παραμέτρου σε ζώνες ίσης τιμής. Η υδραυλική αγωγιμότητα στη x-διεύθυνση χωρίστηκε σε 11 ζώνες, ο συντελεστής οριζόντιας ανισοτροπίας σε 6, η ειδική αποθηκευτικότητα και η ειδική απόδοση σε νερό σε 4 ζώνες η καθεμιά. Οι τιμές των 25 υδρογεωλογικών παραμέτρων προσδιορίστηκαν κατά τη ρύθμιση του μοντέλου. Η ρύθμιση έγινε με τη μέθοδο Gauss-Levenberg-Marquardt με περιορισμούς, όπως αυτή εφαρμόζεται με το PEST, και βασίστηκε σε μηνιαία δεδομένα της περιόδου Σεπτεμβρίου 1992 – Απριλίου 1995. Συγκεκριμένα, χρησιμοποιήθηκαν 1173 παρατηρήσεις φορτίου σε 43 γεωτρήσεις της περιοχής μελέτης και 26 παρατηρήσεις παροχής των πηγών Δίου-Καρίτσας. Το μοντέλο αξιολογήθηκε με ποσοτικά και ποιοτικά κριτήρια (στατιστικές παράμετροι και γραφικές παραστάσεις διαφόρων μορφών των αποτελεσμάτων). Από την αξιολόγηση προέκυψε ότι το μοντέλο προσομοιώνει ικανοποιητικά τη λειτουργία του φυσικού συστήματος. Με την ανάλυση ευαισθησίας ελέγχθηκε η επίδραση κάθε μιας από τις 38 παραμέτρους του μοντέλου (25 υδρογεωλογικές παράμετροι και 13 παράμετροι παροχών εμπλουτισμού). Για την επαλήθευση του μοντέλου χρησιμοποιήθηκαν μηνιαία δεδομένα της περιόδου Μαΐου 1995 – Αυγούστου 1995 (118 παρατηρήσεις φορτίου και 4 παρατηρήσεις παροχής). Η αξιολόγηση της επαλήθευσης βασίστηκε στα ίδια στατιστικά κριτήρια και γραφήματα με αυτά της αξιολόγησης των αποτελεσμάτων της ρύθμισης. Με την επαλήθευση επιβεβαιώθηκαν τα συμπεράσματα της ρύθμισης και της ανάλυσης ευαισθησίας. Στη συνέχεια, το μοντέλο εφαρμόστηκε για την πρόβλεψη του φορτίου του υδροφορέα, για τρία έτη (Σεπτέμβριος 1995 – Αύγουστος 1998), με στόχο τον έλεγχο δύο σεναρίων που αφορούν τις αντλήσεις για την κάλυψη των αρδευτικών αναγκών της περιοχής. Κατά την ανάπτυξη του μοντέλου της Πιερίας, για την επιλογή της κατάλληλης χρονικής διακριτοποίησης, διερευνήθηκε η χρήση σταθερών χρονικών βημάτων και χρονικών βημάτων που σχηματίζουν γεωμετρική πρόοδο, όσον αφορά τον απαιτούμενο υπολογιστικό χρόνο και την εμφάνιση αποκλίσεων στα υπολογιζόμενα φορτία. Από τη διερεύνηση αυτή προέκυψε ότι η χρήση γεωμετρικής προόδου χρονικών βημάτων (εφόσον δεν επηρεάζει την ευστάθεια ή τη σύγκλιση της λύσης) μειώνει δραστικά τον απαιτούμενο χρόνο, παρουσιάζονται όμως αποκλίσεις των υπολογιζόμενων φορτίων από τα φορτία που υπολογίζονται με σταθερό χρονικό βήμα. Η τελική επιλογή εξαρτάται από τις διαστάσεις του προβλήματος, θα πρέπει όμως σε κάθε περίπτωση να ελέγχεται το μέγεθος και η επίδραση αυτών των αποκλίσεων. Επίσης, όπως φάνηκε και από τη σχετική διερεύνηση, στα μοντέλα μη μόνιμης ροής θα πρέπει να ελέγχεται η διάρκεια επίδρασης της αρχικής συνθήκης, όταν αυτή προκύπτει με παρεμβολή σε υπάρχουσες μετρήσεις φορτίου. Είναι πιθανό τα σφάλματα που εισάγονται από τη χρήση μιας τέτοιας αυθαίρετης αρχικής κατανομής του φορτίου να επηρεάζουν τα υπολογιζόμενα φορτία μέχρι το τέλος της προσομοίωσης. Η επίδραση της οριακής συνθήκης γνωστού οριακού φορτίου στο υδατικό ισοζύγιο του υδροφορέα πρέπει να ελέγχεται, ιδιαίτερα κατά τη διάρκεια των εκτελέσεων ρύθμισης, αφού αυτή η οριακή συνθήκη αποτελεί μια ανεξάντλητη πηγή νερού για τον υδροφορέα. Η ανάλυση ευαισθησίας είναι απαραίτητη για να προσδιοριστεί το μέγεθος αλλά και η χωρική κατανομή αυτής της επίδρασης. Η ιδιαιτερότητα αυτού του τύπου οριακής συνθήκης πρέπει να λαμβάνεται υπόψη και στις εκτελέσεις πρόβλεψης του φορτίου. Οι γραφικές παραστάσεις της χρονικής μεταβολής της υπολογιζόμενης εισροής/εκροής από κάθε κελί που ανήκει στο όριο γνωστού φορτίου, μπορούν να βοηθήσουν στην εκτίμηση του οριακού φορτίου στις εκτελέσεις πρόβλεψης. Με βάση τα αποτελέσματα της διερεύνησης για την επίδραση του γνωστού οριακού φορτίου στα υπολογιζόμενα φορτία και στο υδατικό ισοζύγιο του μοντέλου, προτείνεται η ενσωμάτωση της οριακής συνθήκης γνωστού φορτίου στη διαδικασία της ρύθμισης, εφόσον κατά μήκος αυτού του ορίου παρατηρείται το φαινόμενο του αρτεσιανισμού. Στο μοντέλο της Πιερίας αυτό επιτεύχθηκε με τη βοήθεια τεχνητών παρατηρήσεων μηδενικής εισροής από τα κελιά που ανήκουν στο όριο γνωστού φορτίου. Η συμμετοχή αυτής της οριακής συνθήκης στη ρύθμιση του μοντέλου συνέβαλε στην καλύτερη δυνατή προσέγγιση του τρόπου λειτουργίας του φυσικού συστήματος. Επίσης, η χρήση των τεχνητών παρατηρήσεων συνέβαλε και στον ανεξάρτητο προσδιορισμό των παραμέτρων που οι τιμές τους προσδιορίζονται κατά τη ρύθμιση του μοντέλου. Όπως φάνηκε από τις σχετικές διερευνήσεις, ιδιαίτερη προσοχή πρέπει να δίνεται στην επιλογή του κριτηρίου σύγκλισης της μεθόδου επίλυσης της εξίσωσης της υπόγειας ροής, που αφορά το φορτίο. Η βέλτιστη τιμή του κριτηρίου αυτού διαφοροποιείται ανάλογα με τη μέθοδο επίλυσης που επιλέγεται. Για να αποφασιστεί αν η τιμή του κριτηρίου σύγκλισης του φορτίου είναι ικανοποιητική, θα πρέπει πάντα να ελέγχεται και το υδατικό ισοζύγιο του μοντέλου, το σφάλμα του οποίου δεν πρέπει να ξεπερνά το 0.1%. Η συμμετοχή διαφορετικών ειδών παρατηρήσεων στη διαδικασία της ρύθμισης του μοντέλου, βοηθά στον περιορισμό της συσχέτισης μεταξύ των παραμέτρων των οποίων οι τιμές προσδιορίζονται κατά τη ρύθμιση και στην εύρεση μιας μοναδικής λύσης στο αντίστροφο πρόβλημα. Το πλήθος και η ποιότητα των διαθέσιμων δεδομένων για την περιοχή μελέτης είναι καθοριστικής σημασίας για την κατάρτιση ενός μαθηματικού μοντέλου υπόγειων νερών και την αποτελεσματική προσομοίωση της λειτουργίας του φυσικού συστήματος. Για τη ρύθμιση ενός μοντέλου μη μόνιμης ροής, πρέπει να χρησιμοποιηθεί αρκετά μεγάλος αριθμός παρατηρήσεων, με ικανοποιητική για τις διαστάσεις του προβλήματος χωρική και χρονική κατανομή, που να αντιπροσωπεύουν την αντίδραση του συστήματος κάτω από διαφορετικές συνθήκες. Έτσι αυξάνεται η αξιοπιστία των προβλέψεων του μοντέλου. Η αξιολόγηση ενός αυτόματα ρυθμισμένου μοντέλου πρέπει να στηρίζεται σε στατιστικά κριτήρια που αφορούν συνολικά το μοντέλο, τα υπόλοιπα μεταξύ των παρατηρήσεων και των αντίστοιχων υπολογισμένων τιμών, τις απόλυτες τιμές αυτών των υπολοίπων και τις τελικές τιμές των παραμέτρων που προσδιορίστηκαν με τη ρύθμιση. Επίσης απαραίτητη είναι και η παρουσίαση του υδατικού ισοζυγίου του μοντέλου. Η χρήση πολλών κριτηρίων θεωρείται απαραίτητη, γιατί κάθε ένα από αυτά ελέγχει το μοντέλο από διαφορετική σκοπιά και επίσης γιατί τα περισσότερα από αυτά τα κριτήρια προέρχονται από την ανάλυση χρονοσειρών και η εξαγωγή τους βασίζεται σε μεγαλύτερο πλήθος δεδομένων από όσα είναι συνήθως διαθέσιμα στα μοντέλα υπόγειων νερών. Όπως φάνηκε από την ανάλυση όλων των σταδίων, κατά την κατάρτιση μοντέλων υπόγειων νερών με τη βοήθεια σύγχρονων πακέτων λογισμικού, όλες οι διαδικασίες που γίνονται αυτόματα (π.χ. παρεμβολές σε μετρήσεις) και όλες οι προεπιλογές του λογισμικού κατάρτισης πρέπει να ελέγχονται και ο χρήστης να γνωρίζει τον τρόπο λειτουργίας τους και το θεωρητικό τους υπόβαθρο. Το ίδιο ισχύει και για την αυτόματη ρύθμιση ενός μοντέλου με τη βοήθεια εμπορικού ή κοινής κτήσης λογισμικού. Ο χρήστης δεν πρέπει να αρκείται στην απλή εφαρμογή του προγράμματος, βασιζόμενος στις προεπιλεγμένες τιμές των παραμέτρων της μεθόδου βελτιστοποίησης. Οι κατάλληλες τιμές αυτών των παραμέτρων, για κάθε συγκεκριμένο μοντέλο, προκύπτουν μετά από διερεύνηση.