Διερεύνηση των εδαφικών υποχωρήσεων από υπεράντληση υδροφορέων στο δυτικό θεσσαλικό κάμπο, με βάση τη γεωτεχνική συμπεριφορά των σχηματισμών και τη συμβολή γεωστατιστικού μοντέλου

Abstract

Subsidence phenomena due to overexploitation of aquifers have been observed in many regions around the world, but also in the Greek territory. In the latter, these phenomena have been also manifested in the Western plain of Thessaly basin during the last few decades. The occurrence of land subsidence depends on the nature of the formations of the subsoil and particularly of their geotechnical properties. In order to investigate the subsidence mechanism at the above plain, the data regarding lithology and geomechanical properties of soil horizons were examined and analysed, by means of geostatistical methods.Thessaly basin is located in Central Greece and consists a lowland area that is mainly drained by Pinios River. This basin is subdivided by a group of hills in two well distinguished sub-basins, the Western and the Eastern ones. They are two individual hydrologic basins, developing high potential aquifers. According to geological studies of Western Thessaly basin, Mesozoic alpine formations and Molassic formations outcrop in the margins of the study area, while post alpine deposits occupy in the lowland of the basin. The Mesozoic alpine formations consist of schist-cherts, ophiolites, limestones and flysch (Mariolakos et al., 2001; Rozos & Tzitziras, 2002), while the Molassic formations mainly of conglomerates and locally sandstones. The Alpine formations belong to the Pelagonian - Subpelagonian geotectonic zones and they constitute the bedrock of post alpine deposits of the West Thessaly basin.The alluvial deposits of this basin (the Western one) which is the study area of this Thesis constitute a system of unconfined shallow aquifers, extending in the upper layers and successive confined - artesian aquifers developing in the deeper permeable layers (Marinos et al. 1995, 1997). This system apart from the percolated surface water is also supplied by water through the lateral infiltration from the karstic aquifers of the alpine carbonate formations, outcropping in the margins of the basin. In general, the richest aquifers are developed in the western sub-basin of Thessaly plain, due to their rich supply both from the big infiltrating part of surface runoff and lateral infiltration (Kallergis, 1971, 1973; Andronopoulos et al., 1991).The overexploitation of these aquifers during the last decades, followed by a continuous lowering of the ground water level year by year, led to the manifestation of land subsidence phenomena, with extended damages in certain sites. In the Eastern Thessaly basin the land subsidence phenomena occurred for the first time in 1996 (Soulios 1997; Soulios 1980). On the contrary, in Western Thessaly basin (the study area) these phenomena were firstly manifested in 2002 as surface ruptures in the eastern part of it, namely in Farsala and Stavros sites (Rozos et al., 2010). Farsala city is built partly on the bedrock alpine formations and on the Pleistocene deposits. The variations of the geotechnical behaviour for the foundation formations led to the manifestation of numerous tensile fractures, in several sections of the town. Precisely, in the center of the town, an area extending 50 m x 360 m was damaged. For example, the road pavements present multiple fractures, redisplayed afterwards any repair works. Also, several buildings intersected by ruptures were intensively damaged requiring expensive reconstruction works. Small ground ruptures have also been occurred in the northern part of the town and in an area covering 180 m × 200 m. In addition, beyond the south western limits of the town and at the west of the railway line, two more extensive ruptures were observed with total length 1,000 m and 2,500 m respectively. In Stavros small town, the main ground rupture was found westwards the railway line. This tensile rupture has a total length of about 2,100 m, and a vertical displacement at a rate of 60 cm. The trace of the rapture affects road pavements and numerous buildings (Apostolidis & Georgiou, 2007). The buildings founded along the trace of the ruptures present several damages, such as cracks in the stonework, distortions in doors, windows, stockyards and pavements. Furthermore, several ground ruptures are located at the south of the town, intersecting cultivated areas.The overexploitation of the ground water resulted in the activation of the subsidence mechanism of the discharged aquifers and subsequently led to the manifestation of the accompanying phenomena on the surface, apart from the land depression. Therefore, along the margins of the basin where the bedrock outcrops and consequently the thickness of the Quaternary formations are small, ruptures of the ground occur, as a result of the tensile forces action (Rozos et al., 2010). On the contrary, in the parts of the basin with thick Quaternary formations, the compaction of the formations can become noticeable at the first stages mainly by the rising of well pipes from the ground.In order to investigate the mechanism of land subsidence and its future manifestation in some other parts of the basin, data regarding the lithology and the geomechanical characteristics of the involved soil horizons were examined and analysed. The evaluated data, were collected both from several hundreds of sampling boreholes, which had been drilled by local authorities and private companies for various technical works, but also from laboratory tests carried out in the frame of this study (Sideri et al., 2014). Based on the above data but also on an extensive field work, a study of the distinction of the soil horizons regarding their physical characteristics and mechanical properties as well as their classification was conducted. This led in the establishment of a clear relationship between some mechanical properties related to subsidence manifestation and the classification of the distinguished soil horizons.The research results in combination with field work, led to the compilation of two maps. The first one includes the location of all geotechnical and hydrogeological boreholes used in order to assess and evaluate the physical characteristics and mechanical properties of the main engineering geological units. The extended field work along with the evaluated data from these boreholes led to compilation of the second map i.e. the Engineering geological map of the Quaternary deposits of Western Thessaly plain. This map shows the surface development of the eight (8) main lithological - engineering geological units on a topographic background with contour intervals of 200 m especially in the surroundings highlands of the study plain. This map was designed to facilitate the study of land subsidence in the Western Thessaly plain and it is a very useful tool for local public services, as well as individuals involved in urban planning and land development.Apart from the study of the geotechnical characteristics for the separated units of the Quaternary formations and their special development, techniques of multivariate statistics and geostatistics were combined and compared to evaluate the estimation methods of the main mechanical properties, with special reference to compression index (Cc) from oedometer tests. Principal Component Analysis (PCA) was applied in order to highlight the relationships between the geotechnical parameters (Sideri et al., 2014). Through cross-validation analysis, kriging and cokriging were tested as estimators of the Cc.The dataset included geotechnical information of 4,671 samples taken from 813 boreholes. The spatial distribution of the samples with the elevation was uniform, namely 2,338 samples (approximately 55%) between 100 m and 120 m a.s.l. and just 382 samples (approximately 9%) above 150 m a.s.l. The selected samples had a set of geotechnical information comprising the most relevant 21 physical properties - gravel, sand, silt, clay and silty clay fragment percentage %, liquid limit (wl), plastic limit (wp), plasticity index (Ip), clay activity (A), natural water content (wc), consistency index (Ic), bulk unit weight (γb), dry unit weight (γd), particle specific gravity (Gs), void ratio (e), porosity (n), organic matter percentage (%) - and 12 mechanical parameters - max unconfined compressive strength (qu) and corresponding strain (ε) from unconfined compression test, cohesion (c) and friction angle (φ) from shear test, cohesion (c) and friction angle (φ) from triaxial test (UU, CUPP) compression index (Cc) and void ratio from oedometer test (e0). The analysis of the geotechnical information revealed the relationships between soil classification in various units and their mechanical properties. According to USCS, the alluvial deposits in the Western Thessaly plain were subdivided into 28 classes. Those classes were grouped into six (6) engineering geological units, the spatial development of which helped in the examination of the mechanism of subsidence manifestation.The scatter plots of the main geotechnical characteristics (compression index – (Cc) and void ratio (e0)) versus depth as a function of the lithologic/textural classes did not highlight any clear correlation between the classification characteristics and the selected properties of the samples. The lack of a clear relationship between soil classification and geotechnical properties is mainly due to the different scale of observation used for the lithologic description, which is done macroscopically and also due to the geotechnical parameterisation, which is done experimentally (Raspa et al., 2008). Despite that the lithologic description of the various horizons is continuous in a sampling borehole, the punctual character of the samples is restricted in a small area around the borehole and therefore the information gained referred to this small area. This fact is generally eliminated through modelling geological and geotechnical data.Therefore, the correlations did not result to a clear relationship between the examined properties and the soil classification of a sample. Thus, because of the numerous variables involved, the spatial extent of the soil horizons and the complex relationships existing among the physical/mechanical properties are more effectively investigated with geostatistics.Interpolation with cokriging was performed in a 100 km x 75 km x 250 m horizontal orthogonal parallelepiped grid, divided to 1,560,000 elements of 500 m x 500 m x 5 m. 271 samples with Cc value were employed. An isotropic structure with zero nugget effect was used as a model for the variogram. Cokriging was realized by using Cc as target variable with w, e0 and n as auxiliary variables, recorded on a set of 3,021 samples.A corregionalization model was adjusted on the experimental variograms, which were calculated with a lag size of 1 m on the vertical direction and 2,750 m on the horizontal plane. Experimental variograms show an anisotropic variability, which is typical of fluvial sedimentary successions: a small range (8 m) in the vertical direction and a large range (10,100 m) in the horizontal direction. The anisotropy in the range is related to the difference in the sedimentary process along the vertical and horizontal directions. The vertical range of 8 m is also indicative to the dimensions of possible sedimentary bodies, which might have a thickness of about 8 m and much larger horizontal dimensions. The variogram model fitted with the experimental variograms follows the linear corregionalization model. Sill matrixes were automatically adjusted by the software through the conditional fitting of the experimental variograms.The negative Cc values observed in the interpolation results could be explained due to certain inconsistencies between experimental and model covariogram. These inconsistencies cannot be avoided, due to the limitations imposed by the application of positive definiteness to the model. This is the reason why employment of more auxiliary variables for the application of cokriging, could probably result to a very complicated model and is not feasible in general.The complete set of geotechnical information was available for just 219 samples. These samples were analysed to characterize the seven (7) most relevant physical properties and two (2) mechanical parameters. Principal Component Analysis was the first step for subsequent geostatistical analysis. PCA was indeed applied both i) to highlight the relationships among the geotechnical parameters and ii) to synthesize all the information available for kriging application. PCA results showed the sum of the projections of the nine (9) variables on the direction of each of the nine (9) Factors. Due to the redundancy of information related to correlations existing among the parameters, the first two Factors of the PCA explain about 80% of the total variability. Factor 1, which explains 66.91% of the total variability, is highly correlated to Cc, wc, e0, n, γb, γd. Factor 2, which explains 13.02% of the total variability, is correlated to wl, wp and No200. Kriging was applied to the first two Factors, on the same grid as in the previous paragraphs. After that, using the inverse transformation on these two Factors, the Cc value was calculated for each grid element. The cross validation of the results showed a high correlation between true and estimated values of Cc. In addition, unlike cokriging this method does not produce negative Cc values. This is true even though PCA uses eight (8) auxiliary variables for the estimation of Cc, instead of three (3) used in cokriging. The improvement on the model efficiency could be explained as a result of the method’s independence on complicated covariogram structures. The results show that spatial correlation of the target variable is less important than its correlation to the auxiliary variables. Cross-validation demonstrates that kriging with two Principal Component Analysis (PCA) Factors and application of inverse transform gives the best results for Cc estimates (Sideri et al., 2014).Integrating properties, such as relative positions and proportions of different lithofacies, is of highest importance in order to render realistic geological patterns. Plurigaussian Simulation (PS) is an alternative method for conceptual and deterministic modelling for the characterization of lithofacies distribution (Modis & Sideri, 2013). The spatial differentiation of lithofacies was studied in the alluvial aquifer system of West Thessaly basin. For this, PS technique was applied to an extensive set of borehole data from the basin. The main characteristic of the Plurigaussian Simulation technique is that it allows incorporating geological concept in the stochastic simulations. This unique feature of the Plurigaussian Simulation Method is difficult to establish in other similar simulation techniques. The geological rule is not only derived from a statistical analysis of the borehole data but it also allows the inclusion of other forms of geological and empirical knowledge.The PS technique provides a valid indicator cross covariance model because the indicator covariances are directly deduced from the covariance model of the Gaussian variables. Assigning a standard covariance model to these Gaussian variables, guarantees that the indicator covariances are valid covariance functions. The two Gaussian variables are assumed to have stationary covariances. To perform a conditional simulation of the Gaussian variables, experimental multi-Gaussian values are required at the data points. To obtain these values that respect the Gaussian covariance model, a Gibbs sampler algorithm is used. The Gibbs sampler is an iterative algorithm based on Markov chain Monte Carlo (McMC) simulation techniques. The Gaussian values that are assigned to the drill holes are used as conditioning data for the conditional simulation. The two alluvial aquifer systems existing in the Western Thessaly sub-basin can be approximately defined by three boundary surfaces. A thin impermeable layer separates the two aquifer systems, and the lower surface is the limit between the alluvial deposits and the bedrock. The upper aquifer lies between elevations 70 m and 170 m, while the bottom surface of the lower aquifer reaches in depth up to -450 m. The total working area is discretized into 180 x 140 x 71 (1,789,200) blocks of 500 m x 500 m x 10 m. This discretization allows each model block to approximately contain at least one sample. The coordinates of the center of the lower left block are x = 290,000 m, y = 4,330,000 m, z = -450 m.Before starting the stochastic modelling, it is important to take into account the main directions of continuity of the lithology of each horizon because they generally alter in both directions horizontally and vertically (they are not constant in space). Apart from this special arrangement, the base of each horizon imposes a certain kind of spatial continuity. In order to follow this continuity in the simulations, each horizon is deformed to a new shape, with a plane base this time. This deformation is performed by coordinate transformations, while the base of each horizon is used as a reference surface.The lithofacies proportions at each domain must be modelled as a function of the point considered in the field. This modelling procedure is guided mainly by the experimental data. Vertical Proportion Curves (VPCs), first proposed by Matheron et al. (1987), are a simple tool for quantifying the evolution in the amount of each facies or lithotype present as a function of depth. They were computed along lines that are vertical to the chosen reference level.As a first-pass approach, a simple moving average of the experimental indicator values over the study area was considered, in order to adjust the data distribution to the grid. The individual VPCs were used as constraining information to calculate the proportions on the whole 3D grid through kriging estimation. The differentiation between these domains was apparent. A higher proportion of peat-silty clay with organic matter was observable in the second VPC. This was an indication of increased compressibility of the lower aquifer related to the upper one. Similarly, the proportions of sand with gravels were greater in the lower aquifer. Considering the hydraulic conductivity values of the lithofacies in the database, they suggest that the sand with gravels is the most permeable structure. In the PS method, two Gaussian functions G1 and G2 were used. G1 represented the transition between clay - silt grey facies and peat - sand with gravels facies, sand with gravels facies and sand facies, sand facies and basement formations. G2 represented the transition between clay and silt brown facies and all the other hydrofacies.The next step of the PS technique was the inference of the variogram models for the underlying multi-Gaussian functions. Direct adjustment to the experimental variograms was not possible since the only available experimental variograms were the variograms of the indicator functions describing the hydrofacies (one per hydrofacies, plus all the bivariate combinations), while the two variograms needed for the model were the variograms of the underlying and continuous multi-Gaussian functions. The links between all these variograms are complex functions of the truncation process and of the conditioning to the hydrofacies proportions. Therefore, the variogram inference was based on an inverse trial and error procedure in which the ranges of the variograms of the multi-Gaussian fields were adjusted iteratively as follows. Firstly, the type and parameters of the initial variogram models were defined and then the variograms were used to construct an unconditional PS. The variograms of the indicators of the facies were computed numerically from the simulated field. Finally, the initial ranges of the variograms of the multi-Gaussian fields were adjusted manually until an acceptable match was obtained between the experimental and the computed variograms.In the case of West Thessaly basin, the Plurigaussian simulation technique is shown to be effective in reproducing the spatial characteristics of the different lithofacies and their distributions across the study area as determined by the proportions of the borehole samples. Histograms of model vs. experimental lithofacies proportions and indicative cross sections were plotted in order to validate the results. The PS technique was shown to be effective in reproducing the spatial characteristics of the different lithofacies and their distribution across the study area (Modis & Sideri, 2013).In order to investigate the mechanism of land subsidence phenomena and of its future occurrences in the study area, data regarding the temporal variation of ground deformations and water level were examined and analysed. Commonly used methods to measure land subsidence include conventional surveying, borehole extensometers, Global Positioning Systems (GPS), conventional InSAR Interferometric synthetic aperture radar (InSAR) and its recent development PSI technique (Galloway et al, 1999; Maliva & Missimer, 2012). The methods vary in the resolution, cost, and practical spatial density of measurements. Borehole extensometers have the finest resolution and greatest accuracy, but are expensive to be constructed. As a result, there is typically a low spatial density of extensometers in most study areas. Ideally, multiple methods should be used in which, for example, extensometer data may served, for example, as a control regarding less accurate methods with coarser resolutions. In any case, exclusive application of PSI technique in hydrogeological research and monitoring within the past decade, has led to improved results (Ferretti et al. 2000, 2001).Through out of Western Thessaly basin, several water wells were drilled, in order to provide data for various hydro-geological studies (Kallergis 1973; Sogreah S.A. 1974). PSI measurements are also available from the European Space Agency (ESA). The dataset included two types of information: (a) Monthly piezometric time series taken from of 82 water wells from 1974 to 2010 and (b) PSI deformation time series were taken using 68,767 points, from 1992 to 2003. The distribution of water level and ground deformation measurements was not uniform in space and time. The PSI has been applied primarily in urban environments, where the density of stable scatterers typically is quite high (as many as a few hundred per square kilometre). Over natural terrain, the scarcity of stable targets severely limits PSI’s successful application and results in a strong differentiation in sampling density (Galloway & Burbey, 2011).Due to the gaps in the sampling grid as explained above, it was important in the next steps of this work to examine the use of one variable as auxiliary for the estimation of the other. A necessary condition for this is the existence of a certain degree of correlation between the two variables. An initial quantitative examination of the correlation between the piezometric level and the deformation pattern in the PSI data was performed, by linear regression at selected experimental locations. The results did not indicate a clear relation between the above variables. Moderate correlation coefficients between the PSI dataset and the water level fluctuation have been observed, while, in some points, poor correlation coefficients were found.In order to investigate analytically the correlation between deformation and piezometric level measurements in space and time, it was necessary to study the individual variograms and the cross variogram of the above variables. The experimental variograms of ground displacements and water level and also the cross variograms between them were calculated in space and time with a lag size of 1 month on the temporal direction and 1000 m in the spatial domain. These variograms showed an anisotropic variability, which was expected due to the differentiation between the spatial and temporal directions. The sills of the cross variograms were indicative of the correlation between ground deformations and water level in space and time. The above result was in accordance to the theoretically expected: It is well known (Galloway & Burbey, 2011) that there is a direct relation between water level decline and land subsidence. The theoretical variogram fitted to the experimental data followed the linear corregionalization model. Three nested structures, two gaussian and one exponential cosine (Chilès & Delfiner, 2012) were linearly combined to express the spatiotemporal inter-dependencies of ground deformation and water level, individually and between them. The first gaussian model was adjusted to reproduce mainly the spatial continuity of ground deformations and water level individually. This continuity extended isotropically up to 40,000 m. The second gaussian model affects mostly the temporal continuity of both variables, which extended practically up to 10 years. The exponential cosine model with 12 month period simulated the yearly variation of water level due to seasonal fluctuation of rainfall.Since the two variables of interest were correlated, as indicated by their cross variogram, the application of cokriging was expected to be useful for the improvement of estimation. In this direction, both kriging and cokriging were compared for the estimation of ground deformations, with water level as an auxiliary variable. The numerical model was calculated in a spatiotemporal grid (Appiah et al., 2011), extending on 180 km x 140 km in space and 500 months in time, divided to 12,600,000 elements of 500 m x 500 m x 1 month. The coordinates of the center of the lower left block were x= 290,000 m, y= 4,330,000 m and z= 0 which corresponds to December 1971.Study of the variograms of piezometric level and ground deformations in the Western Thessaly basin has shown that it is possible to perform spatiotemporal estimation of these variables, using geostatistics. The cross variogram model indicates that these variables were spatiotemporally correlated. The use of the uniformly sampled water level as an auxiliary variable for the estimation of ground deformations produced more realistic results, by considerably reducing the estimation error (Modis & Sideri, 2015). Model parameters indicated that is safe to interpolate piezometric level within a spatial range of 40 km. It was also possible to predict future values of water level fluctuation. Concerning ground deformations, the variogram model indicated a range of safe interpolation approximately at 40 km. Prediction of future values of ground deformations was also allowable up to 10 years.In order to decrease the estimation errors when forecasting distant values of the model, additional use of other parameters such as soil quality, hydraulic conductivities or various boundary conditions might be of use (Modis & Sideri, 2015). A further future step of the research will be to take into account the above parameters via the incorporation of deterministic laws of ground water flow to the stochastic model.

Το φαινόμενο των εδαφικών υποχωρήσεων εξαιτίας της υπεράντλησης των υπόγειων υδροφορέων, έχει παρατηρηθεί σε πολλές περιοχές ανά τον κόσμο, αλλά και στον Ελληνικό χώρο. Στις περιπτώσεις αυτές, η εκδήλωση ή μη των εδαφικών υποχωρήσεων, εξαρτάται άμεσα από τη φύση των σχηματισμών του υπεδάφους και ειδικότερα από τις γεωτεχνικές τους ιδιότητες. Μία από τις περιοχές στην Ελλάδα στην οποία κατά τις τελευταίες δεκαετίες, το εν λόγω φαινόμενο εμφανίζεται με σοβαρές συνέπειες για τις καλλιέργειες και τις οικιστικές ζώνες είναι και ο Δυτικός Θεσσαλικός κάμπος. Στα πλαίσια της παρούσας διατριβής και προκειμένου να διερευνηθεί η χωροχρονική εκδήλωση του ανωτέρω φαινομένου στην υπόψη περιοχή, εξετάστηκαν και αναλύθηκαν οι λίθοστρωματογραφικές συνθήκες και τα γεωμηχανικά χαρακτηριστικά των εδαφικών οριζόντων, εκτός από τις γνωστές μεθόδους και με τη χρήση της γεωστατιστικής. Αναλυτικότερα, το πρώτο στάδιο της έρευνας περιελάμβανε την αποτύπωση και καταγραφή της παρούσας κατάστασης στον Δυτικό Θεσσαλικό κάμπο. Για τον σκοπό αυτόν, έγινε εκτεταμένη έρευνα πεδίου κατά την οποία αποτυπώθηκαν οι εδαφικές διαρρήξεις και διενεργήθηκαν σταθμημετρήσεις, στο δίκτυο υδρογεωτρήσεων της Περιφέρειας Θεσσαλίας. Στη συνέχεια, πραγματοποιήθηκε χαρτογράφηση και σύνταξη του τεχνικογεωλογικού χάρτη της περιοχής έρευνας. Στον εν λόγω χάρτη, οι τεταρτογενείς αποθέσεις του Δυτικού Θεσσαλικού κάμπου διακρίθηκαν σε οκτώ (8) τεχνικογεωλογικές ενότητες. Μετά την εξαγωγή των πρώτων συμπερασμάτων, τα δεδομένα των γεωτεχνικών γεωτρήσεων και υδρογεωτρήσεων (περισσότερες από 1.000) που έχουν ανορυχθεί στην περιοχή έρευνας, χρησιμοποιήθηκαν για την αποτύπωση της λιθολογικής σύστασης και δομής, καθώς και την προσέγγιση του φαινομένου των εδαφικών υποχωρήσεων, με μεθόδους γεωστατιστικής. Η χρήση της γεωστατιστικής προκύπτει από την ανάγκη αποδοτικότερης χρήσης των δεδομένων, με στόχο την ελαχιστοποίηση του σφάλματος, σε συνθήκες υψηλής αβεβαιότητας. Εξάλλου, η δομή των τεταρτογενών αποθέσεων της λεκάνης του Δυτικού Θεσσαλικού κάμπου χαρακτηρίζεται από έντονη πολυπλοκότητα η οποία, σε συνάρτηση και με την έκτασή της, καθιστά δύσκολη την πλήρη ανάλυση αυτής μόνο με τα υπάρχοντα στοιχεία. Συνεπώς, στο δεύτερο στάδιο της έρευνας, κρίθηκε απαραίτητη η χρήση γεωστατιστικών μεθόδων για την απεικόνιση της τρισδιάστατης λιθολογικής σύστασης και δομής, με σκοπό τον εντοπισμό νέων περιοχών στις οποίες αναμένεται υπό ορισμένες συνθήκες να εκδηλωθούν φαινόμενα εδαφικών υποχωρήσεων. Για τη χωρική χαρτογράφηση των γεωτεχνικών παραμέτρων του εδάφους χρησιμοποιήθηκαν δεδομένα από τις γεωτεχνικές γεωτρήσεις, οι οποίες έχουν ανορυχθεί στην περιοχή. Ειδικότερα, από τα γεωμηχανικά χαρακτηριστικά των σχηματισμών διερευνήθηκε η χωρική μεταβολή του δείκτη συμπιεστότητας (Cc), αυξημένες τιμές του οποίου σε συγκεκριμένους εδαφικούς ορίζοντες, μπορεί να προκαλέσουν τη συμπίεση αυτών και κατά συνέπεια την υποχώρηση του συνόλου των εδαφικών σχηματισμών. Η χωρική μεταβολή του Cc εξετάστηκε σε σχέση με άλλες γεωτεχνικές παραμέτρους των τεταρτογενών αποθέσεων, με την εφαρμογή των μεθόδων cokriging και ανάλυσης κύριων συνιστωσών. Η συσχέτιση φυσικών χαρακτηριστικών και μηχανικών ιδιοτήτων στη λεκάνη του Δυτικού Θεσσαλικού κάμπου με την εφαρμογή των ανωτέρω μεθόδων, οδήγησε πράγματι στον εντοπισμό περιοχών στις οποίες επικρατούν οι συμπιεστοί σχηματισμοί. Ωστόσο, με την απουσία γεωτεχνικών δεδομένων πέραν του βάθους των 30 - 45 m (βάθος των γεωτεχνικών γεωτρήσεων που χρησιμοποιήθηκαν), δεν ήταν δυνατός ο υπολογισμός της τιμής του δείκτη συμπιεστότητας στους βαθύτερους εδαφικούς ορίζοντες της περιοχής έρευνας. Για την αντιμετώπιση του προβλήματος αυτού, αποδόθηκε τιμή συμπιεστότητας σε κάθε λιθολογική ενότητα των βαθύτερων εδαφικών σχηματισμών, η οποία παρουσιάζει παρόμοια γεωτεχνική συμπεριφορά με τους ορίζοντες που διακρίθηκαν έως το βάθος διερεύνησης. Ως εκ τούτου, για την ορθότερη απεικόνιση της τρισδιάστατης λιθολογικής σύστασης και δομής των τεταρτογενών αποθέσεων, αλλά και την αποφυγή εξομαλύνσεων, εφαρμόστηκε η μέθοδος της πολλαπλά Γκαουσιανής προσομοίωσης. Από τα αποτελέσματα της προσομοίωσης επιβεβαιώνεται ότι η λιθολογική σύσταση και η δομή, η οποία εκφράζεται και από την κοκκομετρική διαβάθμιση των εδαφικών σχηματισμών στην εν λόγω περιοχή, παρουσιάζουν διαφοροποίηση στο ΒΔ/κό συγκριτικά με το ΝΑ/κό τμήμα της λεκάνης του Δυτικού Θεσσαλικού κάμπου. Πιο συγκεκριμένα, το μοντέλο επιβεβαιώνει την ύπαρξη μεγαλύτερου ποσοστού αδρομερών σχηματισμών στο ΒΔ/κό συγκριτικά με το ΝΑ/κό τμήμα αυτής. Η συγκεκριμένη διαφοροποίηση επαληθεύεται από έρευνες οι οποίες έχουν πραγματοποιηθεί κατά το παρελθόν στην εν λόγω περιοχή. Μετά τον εντοπισμό των περιοχών οι οποίες δύνανται υπό συνθήκες να υποστούν εδαφική υποχώρηση, διερευνήθηκε η ποσοτική εξάπλωση του φαινομένου στον χώρο και τον χρόνο, σε σχέση με τη γενεσιουργό αιτία, δηλαδή την πτώση στάθμης του υπόγειου υδροφορέα. Για τον σκοπό αυτόν, χρησιμοποιήθηκε χρονοσειρά δορυφορικών μετρήσεων από το ερευνητικό πρόγραμμα Terrafirma του Ευρωπαϊκού Οργανισμού Διαστήματος, βάσει της τεχνικής της Συμβολομετρίας Σταθερών Ανακλαστήρων. Από τα αποτελέσματα των δορυφορικών μετρήσεων εντοπίστηκαν οι περιοχές οι οποίες παρουσιάζουν έντονες παραμορφώσεις. Με βάση τα δεδομένα των δορυφορικών μετρήσεων και τη χρήση των μεθόδων kriging και cokriging εκτιμήθηκαν οι αλληλοσχετιζόμενες χωροχρονικές μεταβολές των μετακινήσεων του εδάφους και της στάθμης του υπόγειου νερού και πραγματοποιήθηκε πρόβλεψη των τιμών τους. Συμπερασματικά, η διερεύνηση των μεταβολών της πιεζομετρίας και των μετακινήσεων του εδάφους στον Δυτικό Θεσσαλικό κάμπο κατέδειξε ότι, με την εφαρμογή της γεωστατιστικής είναι δυνατή η χωροχρονική εκτίμηση αυτών εντός συγκεκριμένων ορίων. Η χρήση της πιεζομετρίας ως βοηθητικής μεταβλητής για την εκτίμηση και πρόβλεψη των μετακινήσεων του εδάφους, αποδείχθηκε ότι οδηγεί σε περισσότερο ρεαλιστικά αποτελέσματα, μειώνοντας σημαντικά το σφάλμα εκτίμησης. Τα αποτελέσματα της διατριβής, όσον αφορά τόσο στην απόδοση της επικινδυνότητας, όσο και στη χωροχρονική εκτίμηση του μεγέθους των συνεπειών, αναμένεται να συμβάλλουν στον σχεδιασμό των δράσεων ανάσχεσης της εκδήλωσης των υποχωρήσεων του εδάφους, όπου αυτές έχουν σημειωθεί, αλλά κυρίως στη λήψη προληπτικών μέτρων στις περιοχές της λεκάνης οι οποίες είναι μεν επιρρεπείς σε εδαφική υποχώρηση, πλην όμως οι υπεραντλήσεις δεν έχουν φθάσει στον βαθμό εκείνον που είναι ικανός να οδηγήσει σε ορατές ενδείξεις στην επιφάνεια, όπως εδαφικές διαρρήξεις κ.λπ.