Υδροδυναμική και περιβαλλοντική προσομοίωση του Παγασητικού κόλπου

Abstract

This Phd thesis studies and investigates the hydraulic -hydrodynamic simulation of Pagasitikos gulf using a three-dimensional computational hydrodynamic model in order to provide useful and accurate information and predictions about the behavior of the gulf, which cannot be obtained from simple measurements and monitoring. This model takes into account and simulates in the best possible way data such as the actual precise topography and bathymetry of Pagasitikos gulf, meteorological data for the wider area such as wind speed, wind direction, atmospheric pressure, relative humidity, rainfall, temperature, and tides. For the model construction simulations, Elcom 2.2 (Estuary, Lake and Coastal Ocean Model) software was used, available at the Hydraulics Department of DUTh, plus the later version of the program AEM 3D (Aquatic Ecosystem Model), which was purchased with an annual subscription license from Hydronumerics in Australia. ELCOM uses hydrodynamic and thermodynamic models to simulate the temporal behaviour of stratified flows under the influence of environmental actions.The numerical method applied by ELCOM solves the non-steady Navier-Stokes equations for incompressible flow, ignoring the non-hydrostatic pressure terms. For horizontal turbulence, the Reynolds-averaged Navier-Stokes equations are solved and turbulence closure is performed using a turbulent viscosity coefficient, while for vertical turbulence, a vertical mixing layer model is applied. Furthermore, the transport of passive and active tracers (e.g., tracers, salinity, and temperature) is performed using the conservative semi-selective discretization method ULTIMATE QUICKEST. The ELCOM model has been successfully applied in other related studies in Greece, such as in the north Aegean Sea and the Gulf of Evia, where tides were also simulated. A simulation field was selected and created, enclosed by the coastlines of the Pagasitikos gulf and north Evia with two open sea cross-sections on the left and right, south of the Trikeri channel, where two open boundary conditions are simulated to the east and west, with lengths of 10 km and 5.2 km, respectively. The first open sea section located in the Artemisio Strait (east of the Maliakos gulf) and the second is eastwars at the narrowest point between northern Evia and the southern coastline of Pagasitikos (west of Skiathos). The maps that were digitized are the entire Pagasitikos gulf and a section of northern Evia. This area provides ideal simulation conditions as it has no large open sea boundaries, is enclosed by solid boundaries with two relatively small and defined open sections where boundary conditions are very well applied. Taking into account the morphology of the bay, the size of the surface area, and the final shape of the computational domain (47.600 x 31.600 m) with a surface area of 1.504.160 m2, a fixed computational cell size of (400 x 400 m) was selected. The maximum simulation depth reaches 105 m and 41 horizontal layers were selected from the surface to the deepest point of the bottom with a fixed height of dz = 2.5 m, creating 385,441 computational cells. Because the accuracy and stability of the simulation is negatively affected by the degree of variation in the height of the horizontal layers (this is described as non-uniformity of the layers), a three-dimensional grid with a uniform layer thickness (uniform dz grid) was constructed to maximize the accuracy of the simulation. Data from two sources were used: a) the Copernicus - Mediterranean Sea Waves Analysis and Forecast system, specifically data from the MEDSEA_ANALYSISFORECAST_WAV_006_017 file, which contains measurements of sea level and wave height from May 4, 2019, to the present, and b) the Poseidon system of HELLENIC CENTRE FOR MARINE RESEARCH, which has systematic measurements, among other field data, and values for sea level from the measuring station in Skopelos. For a full calendar year, the sea level difference was entered at the two open boundaries, and the model automatically generated the velocity field. Complete meteorological data from the National Observatory of Athens were used for atmospheric loads. These are measurements from the nearest station installed in Volos, recorded every 10 minutes from 2007 to 2018. The model includes the Coriolis force due to the Earth's rotation. A data storage grid was selected and created at the following depths (layers): surface (0 m), 2.5 m, 5 m, 7.5 m, 10 m, 12.5 m, 15 m, 20 m, 25 m, 35 m, 45 m, 55 m, 75 m, 95 m, and at the bottom. Three vertical layers of computational cells (curtains) were defined to store the simulation results. Two of these levels run west-east, one at the height of Almyros and one at the Trikeri strait, and a third runs north-south through the port of Volos. In addition, 24 points were defined on the coastline at the location of organized beaches. In order to select the optimal time step (dt), a sensitivity analysis was performed using short-duration test simulations for 60 sec, 120 sec, 180 sec, and 240 sec time steps. The average values of the measurements were used and the measured velocities were compared with the velocities from simulations with the above mentioned time steps dt. The comparison of model velocity values for the four time steps at various points led to the selection of the optimal time step dt =180 sec. The validity and reliability of the results of this simulation are confirmed by published results of measurements taken from samples. To our knowledge, there is no other corresponding computational simulation for the area under consideration. This study provides visualisations of a series of hydrodynamic and qualitive quantities at the surface and at various depths and investigates the mechanisms of water renewal in the bay waters. The simulation was carried out for one year and the resulting values for velocity, temperature, density, salinity, and water residence time were depicted on the surface, at various depths and in vertical sections of the bay. The results were analyzed and evaluated, and the way in which bathymetry, topography, tides, and atmospheric loading determine and influence hydrodynamic circulation and the mixing of water masses. The following were examined in detail: a) The speed, temperature, density, and salinity of the water in Pagasitikos gulf during a one-year simulation at the surface and throughout its entire depth. b) The age, movement, residence, and renewal of the waters of Pagasitikos gulf with open sea water in all computational cells. c) The effect of density stratification on the limitation of vertical mixing and the entrapment of marine masses. d) The creation and movement of marine currents in the semi-enclosed sea basin of the study area. e) The communication and exchange of water between Pagasitikos gulf and north Evian gulf. f) The movement of sea currents between north Evia and Pagasitikos gulf through the two open sea cross-sections. g) The inflow and outflow of water from the gulf, in order to obtain a complete picture of the hydrodynamic behaviour of Pagasitikos gulf. The study and investigation of hydrodynamic circulation within Pagasitkosc gulf through computational simulation, as carried out in this doctoral thesis, contributes to research and science in the following ways: 1) It provides a detailed presentation and assessment of the hydrodynamic characteristics of the water column of Pagasitikos gulf at specific critical points and areas. It identifies areas that, due to their slower water renewal compared to others, are more prone to pollution from urban wastewater discharged into the wider bay area. 2) It highlights and identifies the communication and water exchange between Pagasitikos gulf and north Evian gulf. 3) It contributes to the understanding of the behavior of the bay and accurately identifies the three-dimensional circulation of water under prevailing natural conditions and loads. 4) It provides a reliable numerical model, through the confirmation and evaluation of results, which can assist in the broader environmental planning of the area and in preventing environmental degradation of the bay. 5) The proposed numerical investigation methodology can serve as the basis for a comprehensive and reliable research tool for longer-term simulations to predict the future behavior of the bay.6) It enables more focused investigation and simulation of specific conditions. The simulation can focus on smaller parts of the area, such as the mouths of streams and rivers, or simulate extreme events such as floods and others.

Η παρούσα διδακτορική διατριβή μελετά και διερευνά την υδραυλική, υδροδυναμική και κατ’ επέκταση περιβαλλοντική προσομοίωση του Παγασητικού κόλπου με τη χρήση τρισδιάστατου υπολογιστικού υδροδυναμικού μοντέλου με σκοπό να δώσει χρήσιμες και ορθές πληροφορίες και προβλέψεις για την συμπεριφορά του κόλπου, οι οποίες δεν μπορούν να προκύψουν από απλές μετρήσεις και παρακολούθηση. Στο μοντέλο αυτό λαμβάνονται υπόψη και προσομοιώνονται με τον βέλτιστο τρόπο δεδομένα όπως η πραγματική ακριβής τοπογραφία και βυθομετρία του Παγασητικού κόλπου και τα μετεωρολογικά δεδομένα της ευρύτερης περιοχής όπως η ταχύτητα ανέμου, η κατεύθυνση του ανέμου, η ατμοσφαιρική πίεση, η σχετική υγρασία, το ύψος βροχόπτωσης, η θερμοκρασία και η παλίρροια. Για την κατασκευή του μοντέλου και τις προσομοιώσεις αποφασίστηκε να χρησιμοποιηθεί το λογισμικό Elcom 2.2 (Εstuary, Lake and Coastal Ocean Model), που διαθέτει ο Τομέας Υδραυλικής του ΔΠΘ και η μεταγενέστερη επίσης έκδοση του προγράμματος με την ονομασία AEM 3D (Aquatic Ecosystem Model), το οποίο αγοράστηκε με ετήσια συνδρομητική άδεια χρήσης από την εταιρεία Hydronumerics στην Αυστραλία. Το ELCOM χρησιμοποιεί υδροδυναμικά και θερμοδυναμικά μοντέλα για την προσομοίωση της χρονικής συμπεριφοράς στρωματοποιημένων ροών, υπό την επιρροή περιβαλλοντικών δράσεων. Η αριθμητική μέθοδος που εφαρμόζεται από το ELCOM επιλύει τις μη μόνιμες εξισώσεις Navier-Stokes για ασυμπίεστη ροή, αγνοώντας τους μη υδροστατικούς όρους πίεσης. Για την οριζόντια τύρβη λύνονται οι μεσοποιημένες (μέσες χρονικά) κατά Reynolds εξισώσεις Navier-Stokes και το κλείσιμο της τύρβης πραγματοποιείται μέσω ενός συντελεστή τυρβώδους συνεκτικότητας, ενώ όσον αφορά την κατακόρυφη τύρβη εφαρμόζεται ένα μοντέλο στρώματος κατακόρυφης ανάμειξης. Επίσης, η μεταφορά των παθητικών και ενεργητικών βαθμωτών μεγεθών (π.χ. δείκτες - tracers, αλατότητα και θερμοκρασία) πραγματοποιείται χρησιμοποιώντας τη συντηρητική μέθοδο ημιεπιλεγμένης διακριτοποίησης ULTIMATE QUICKEST. Το μοντέλο ELCOM έχει εφαρμοστεί με επιτυχία και σε άλλες σχετικές έρευνες στον Ελλαδικό χώρο, όπως για παράδειγμα στη θάλασσα του βορείου Αιγαίου και στον Ευβοϊκό κόλπο, όπου μάλιστα προσομοιώθηκε και η παλίρροια. Επιλέχθηκε και δημιουργήθηκε ένα πεδίο για προσομοίωση, το οποίο περικλείεται από τις ακτογραμμές του Παγασητικού κόλπου και της βόρειας Εύβοιας, με δύο ανοιχτές θαλάσσιες διατομές αριστερά και δεξιά, νότια του διαύλου Τρικερίου, όπου προσομοιώνονται δύο αντίστοιχες ανοικτές οριακές συνθήκες ανατολικά και δυτικά μήκους 10 km και 5,2 km αντίστοιχα. Αυτές βρίσκονται η πρώτη στο Στενό Αρτεμισίου (ανατολικά του Μαλιακού κόλπου) και η δεύτερη ανατολικά στο στενότερο σημείο μεταξύ βόρειας Εύβοιας και νότιας ακτογραμμής Παγασητικού (δυτικά της Σκιάθου). Τα τμήματα των χαρτών που ψηφιοποιήθηκαν είναι ολόκληρος ο Παγασητικός κόλπος, και ένα τμήμα της βόρειας Εύβοιας. Το πεδίο αυτό παρέχει ιδανικές συνθήκες προσομοίωσης καθώς δεν έχει μεγάλα ανοιχτά θαλάσσια όρια και περικλείεται από στερεά όρια με δύο σχετικά μικρά και οριοθετημένα ανοιχτά τμήματα όπου εφαρμόζονται πολύ καλά οριακές συνθήκες. Λαμβάνοντας υπόψη την μορφολογία του κόλπου, το μέγεθος της επιφάνειας και την τελική μορφή του υπολογιστικού πεδίου διαστάσεων (47.600 x 31.600 m) και επιφανείας 1.504.160 m2, επιλέχθηκε σταθερή διάσταση υπολογιστικών κελιών (400 x 400 m). Το μέγιστο βάθος προσομοίωσης φτάνει τα 105 m και επιλέχθηκαν 41 οριζόντια επίπεδα (layers) από την επιφάνεια έως το βαθύτερο σημείο του πυθμένα σταθερού ύψους dz = 2.5 m και δημιουργήθηκαν 385.441 υπολογιστικά κελιά. Επειδή η ακρίβεια και η σταθερότητα της προσομοίωσης επηρεάζεται αρνητικά από τον βαθμό μεταβολής του ύψους των οριζόντιων επιπέδων (αυτό περιγράφεται ως non-uniformity of the layers), κατασκευάστηκε ένα τρισδιάστατο πλέγμα με σταθερό πάχος στρώσεων (uniform dz grid) ώστε να μεγιστοποιείται η ακρίβεια της προσομοίωσης. Χρησιμοποιήθηκαν δεδομένα από δύο πηγές: α) το σύστημα Copernicus - Mediterranean Sea Waves Analysis and Forecast και συγκεκριμένα δεδομένα από το αρχείο MEDSEA_ANALYSISFORECAST_WAV_006_017, όπου υπάρχουν μετρήσεις για το ύψος της στάθμης θάλασσας και κυμάτων από 4/5/2019 μέχρι σήμερα και β) το σύστημα Poseidon του ΕΛΚΕΘΕ το οποίο έχει συστηματικές μετρήσεις, μεταξύ των άλλων δεδομένων πεδίου, και τιμές για την στάθμη της θάλασσας από τον σταθμό μέτρησης στη Σκόπελο. Για ένα πλήρες ημερολογιακό έτος έγινε εισαγωγή της διαφοράς στάθμης της θάλασσας στα δύο ανοικτά όρια, οπότε το μοντέλο παρήγαγε αυτόματα το πεδίο ταχυτήτων. Για τις ατμοσφαιρικές φορτίσεις χρησιμοποιήθηκαν πλήρη μετεωρολογικά δεδομένα, προερχόμενα από το ΙΕΠΒΑ του Εθνικού Αστεροσκοπείου Αθηνών. Πρόκειται για μετρήσεις από τον πλησιέστερο σταθμό που έχει εγκατασταθεί στον Βόλο, ανά 10 min καταγεγραμμένες από το 2007 μέχρι το 2018. Το μοντέλο συμπεριλαμβάνει την δύναμη Coriolis, λόγω περιστροφής της γης. Επιλέχθηκε και δημιουργήθηκε ένα πλέγμα αποθήκευσης δεδομένων στα εξής βάθη (layers): επιφάνεια (0 m), 2.5 m, 5 m, 7.5 m, 10 m, 12.5 m, 15 m, 20 m, 25 m, 35 m, 45 m, 55 m, 75 m, 95 m και στον βυθό. Για την αποθήκευση αποτελεσμάτων της προσομοίωσης ορίσθηκαν τρία κατακόρυφα επίπεδα υπολογιστικών κελιών (curtains). Τα δύο από αυτά επίπεδα έχουν κατεύθυνση δύση-ανατολή, ένα στο ύψος του Αλμυρού και ένα στο στενό Τρίκερι και ένα τρίτο με κατεύθυνση βορά-νότο διέρχεται μέσα από το λιμάνι του Βόλου. Επιπλέον ορίσθηκαν 24 σημεία στην ακτογραμμή στη θέση που βρίσκονται οργανωμένες παραλίες. Προκειμένου να επιλεγεί το βέλτιστο χρονικό βήμα (dt), έγινε ανάλυση ευαισθησίας με δοκιμαστικές προσομοιώσεις μικρής διάρκειας, για διάφορες τιμές χρονικού βήματος. Χρησιμοποιήθηκαν οι μέσες τιμές των μετρήσεων και έγινε σύγκριση των καταμετρημένων ταχυτήτων με τις ταχύτητες από προσομοιώσεις με dt 60 sec, 120 sec 180 sec και 240 sec. Η σύγκριση των τιμών ταχυτήτων μοντέλου για τα τέσσερα χρονικά βήματα σε διάφορα σημεία, οδήγησε στην επιλογή του βέλτιστου χρονικού βήματος dt = 180 sec. Η εγκυρότητα και η αξιοπιστία των αποτελεσμάτων της παρούσας προσομοίωσης επιβεβαιώνεται από δημοσιευμένα αποτελέσματα μετρήσεων πραγματοποιηθέντων δειγματοληψιών. Από την παρούσα έρευνα προκύπτουν απεικονίσεις μεγεθών στην επιφάνεια και σε διάφορα βάθη και διερευνώνται οι μηχανισμοί ανανέωσης των νερών του κόλπου. Η προσομοίωση έγινε για ένα έτος και οι τιμές ταχύτητας, θερμοκρασίας, πυκνότητας, αλατότητας και χρόνου παραμονής των υδάτων που προέκυψαν, απεικονίστηκαν στην επιφάνεια και σε κατακόρυφες τομές του κόλπου. Έγινε ανάλυση και αξιολόγηση των αποτελεσμάτων και διερευνήθηκε λεπτομερώς ο τρόπος που η βυθομετρία, η τοπογραφία, η παλίρροια και η ατμοσφαιρική φόρτιση καθορίζουν και επηρεάζουν την υδροδυναμική κυκλοφορία καθώς και την ανάμιξη των υδάτινων μαζών. Εξετάστηκαν αναλυτικά τα εξής: α) Η ταχύτητα, η θερμοκρασία, η πυκνότητα, η αλατότητα του νερού του Παγασητικού κόλπου στη διάρκεια ενός έτους προσομοίωσης στην επιφάνεια αλλά και σε όλο το βάθος. β) Η ηλικία, η κίνηση, η παραμονή και ανανέωση των νερών του Παγασητικού κόλπου με το νερό του ανοιχτού πελάγους σε όλα τα υπολογιστικά κελιά. γ) Η επίδραση της στρωμάτωσης πυκνότητας στον περιορισμό της κατακόρυφης ανάμιξης και στον εγκλωβισμό των θαλάσσιων μαζών. δ) Η δημιουργία και κίνηση των θαλάσσιων ρευμάτων στην ημίκλειστη θαλάσσια λεκάνη της περιοχής μελέτης. ε) Η επικοινωνία και ανταλλαγή νερού του Παγασητικού κόλπου με τον βόρειο Ευβοϊκό κόλπο.στ) Η κίνηση των θαλάσσιων ρευμάτων μεταξύ βόρειας Εύβοιας και Παγασητικού ανάμεσα στις δύο ανοιχτές θαλάσσιες διατομές και τέλος. ζ) Η είσοδος και έξοδος του νερού από τον κόλπο, ώστε να υπάρξει μια ολοκληρωμένη εικόνα της υδροδυναμικής συμπεριφοράς του Παγασητικού κόλπου. Η μελέτη και διερεύνηση της υδροδυναμικής κυκλοφορίας στο εσωτερικό του Παγασητικού κόλπου μέσω υπολογιστικής προσομοίωσης, όπως πραγματοποιείται στην παρούσα διδακτορική διατριβή, συμβάλει ερευνητικά και επιστημονικά στα εξής: 1) Παρέχει μια λεπτομερή παρουσίαση και εκτίμηση των υδροδυναμικών χαρακτηριστικών της υδάτινης στήλης του Παγασητικού κόλπου σε συγκεκριμένα κρίσιμα σημεία και περιοχές. Καθορίζονται περιοχές που, λόγω της πιο αργής ανανέωσης των υδάτων τους σε σύγκριση με άλλες, είναι πιο επιρρεπείς στη ρύπανση από τα αστικά λύματα που απορρίπτονται στην ευρύτερη περιοχή του κόλπου. 2) Αναδεικνύει και προσδιορίζει την επικοινωνία και ανταλλαγή νερού μεταξύ Παγασητικού και βόρειου Ευβοϊκού κόλπου. 3) Συμβάλει στην κατανόηση της συμπεριφοράς κόλπου και προσδιορίζει με ακρίβεια την τρισδιάστατη κυκλοφορία του νερού υπό τις επικρατούσες φυσικές συνθήκες και φορτίσεις. 4) Παρέχει ένα αξιόπιστο αριθμητικό μοντέλο, μέσω της επιβεβαίωσης και αξιολόγησης των αποτελεσμάτων, που μπορεί να βοηθήσει στον ευρύτερο περιβαλλοντικό σχεδιασμό της περιοχής και στην πρόληψη της περιβαλλοντικής υποβάθμισης του κόλπου. 5) Η αριθμητική μεθοδολογία διερεύνησης που προτείνεται, μπορεί να αποτελέσει τη βάση αλλά και ένα ολοκληρωμένο και αξιόπιστο εργαλείο έρευνας για μεγαλύτερες χρονικά προσομοιώσεις για πρόβλεψη της μελλοντικής συμπεριφοράς του κόλπου. 6) Δίνει την δυνατότητα πιο εστιασμένης διερεύνησης και προσομοίωσης εξειδικευμένων συνθηκών και φορτίσεων. Μπορεί να εστιαστεί η προσομοίωση σε μικρότερα τμήματα της περιοχής, όπως σημεία εκβολής χειμάρρων και ποταμών ή να προσομοιωθούν ακραία φαινόμενα όπως πλημμύρες και άλλα.