Abstract

In the current context of climate change, improving groundwater monitoring and management is an important issue for human communities in arid environments. The exploitation of groundwater resources can trigger land subsidence producing damage in urban structures and infrastructures. Alto Guadalentín aquifer system in SE Spain has been exploited since 1960 producing an average piezometric level drop of 150 m. This work presents a groundwater model that reproduces groundwater evolution during 52 years with an average error below 10%. The geometry of the model was improved introducing a layer of less permeable and deformable soft soils derived from InSAR deformation and borehole data. The resulting aquifer system history of the piezometric level has been compared with ENVISAT deformation data to calculate a first-order relationship between groundwater changes, soft soil thickness, and surface deformation. This relationship has been validated with the displacement data from ERS and Cosmo-SkyMed satellites. The resulting regression function is then used as an empirical subsidence model to estimate a first approximation of the deformation of the aquifer system since the beginning of the groundwater extraction, reaching 1 to 5.5 m in 52 years. These rough estimations highlight the limitations of the proposed empirical model, requiring the implementation of a coupled hydrogeomechanical model.

1. Introduction

Aquifer overexploitation and groundwater reservoirs depletion represent a problem impacting large areas worldwide [1, 2], especially over arid and semiarid areas. Agricultural and urban use of aquifers resources, especially during dry seasons, has led to water mining which often leads to reservoir depletion and saline intrusion [3, 4]. These problems have taken more relevance in a global climate change context that predicts intensive dry seasons over vulnerable areas [5, 6]. Ground subsidence, accompanying aquifer system depletion, is a common hazard that has important social and economic repercussions [7]. The effect of piezometric level changes over aquifer systems and associated consolidation can be explained by Terzaghi’s effective stress principle [8]. This type of soil deformation may produce building and infrastructure damage [9, 10], changes in flooding risk areas, or loss of water storage capacity [11]. In order to prevent major problems, a good knowledge of groundwater resources changes and effects should be a main priority for water management entities.

The use of numerical groundwater flow models to reproduce the hydrodynamic behavior of an area and simulate spatial and temporal changes in piezometric levels, together with synthetic aperture radar differential interferometry (DInSAR) that detects surface movements, allows a global approach to assess main aquifer threats. Both techniques have been used to understand subsidence related to groundwater withdrawal [12] and calibrate hydrogeological and geomechanical parameters [13].

This study focuses on the Alto Guadalentín basin located in the Mediterranean arc, SE of Spain. The Spanish Mediterranean arc is a drought vulnerable area that experienced three important dry periods between 1990 and 2012 [14], where 68 aquifer systems have been declared partially overexploited and 10 as completely overexploited by national authorities. Agriculture, traditionally the most important economic activity in the area, is being progressively replaced by urban and touristic activities. These activities still have an important impact over groundwater resources.

Piezometric data from the Alto Guadalentín basin show a continuously descending piezometric level over the last 60 years from a few meters below land surface in the 60s to approximately 200 m below land surface nowadays. Previous studies have used multi-sensor SAR images from ERS, ENVISAT, ALOS, and Cosmo-SkyMed (CSK) satellites to obtain a 20-year (1992–2012) deformation map over the area [1517]. The results have revealed one of the highest subsidence rates measured in Europe (>10 cm/yr).

This paper presents a calibrated numerical groundwater flow model generated using MODFLOW [18] that covers the period from 1960 to 2012 and simulates the spatial and temporal evolution of piezometric level. Subsequently, SAR-derived deformation series has been compared with piezometric change data and soft soils thickness data in order to evaluate a first-order relationship between them. This is the first time that results of a calibrated transient groundwater model, simulating more than 50 years of the Alto Guadalentín aquifer behavior, have been analyzed to understand and quantify the associated subsidence. Moreover, subsidence measurements derived from SAR data have contributed to improving the conceptual and physical groundwater model.

2. Study Area

The Alto Guadalentín basin (Figure 1) is an intramontane sedimentary NE-SW-oriented basin limited by Torrecilla range to the North-West and Almenara range to the South-East and contains an aquifer system that is hydraulically connected to the Bajo Guadalentín aquifer system. Altitudes in the basin range from 300 to 400 m a.s.l. while the surrounding mountains range from 500 to 950 m a.s.l. This basin is an alpine orogenic tectonic depression with a Paleozoic preorogenic metamorphic basement that covers an area of approximately 250 km2. Basement shows a horst and graben pattern with maximum depths of 1000 m below land surface [19]. The Guadalentín depression is bounded to the north by the active Alhama de Murcia fault [20]. The aquifer system developed over the metamorphic basement is constituted by sedimentary materials accumulated during the basin formation [21]. Deepest materials are mostly composed by Miocene marls, acting as impervious seal of the aquifer system, with conglomerates and sandstones on the top of the series. Taking into account the horst and graben pattern of the basement, its thickness has a remarkable variability oscillating from 300 to 900 m depending on the location [22, 23]. Plio-quaternary sediments compose the main and most productive layer of the aquifer system (Figure 2). This layer comprises sediments from near ranges alluvial fans generating sand and gravel lenses embedded in a clay and silt matrix. The upper part of the aquifer is unconfined, while deeper areas have a semiconfined behavior.

The Guadalentín River is the only permanent watercourse in the area (Figure 1), slightly crossing the basin. It presents very low flow rates (average flow rate 0.1 m3/s) during long dry periods with extreme flood events with flows as much as 2000 m3/s [24]. Other occasional watercourses are the Nogalte, Torrecilla, and Béjar ramblas, active after intense rains. The Spanish Mediterranean arc is annually affected by Isolated Depressions at High Levels [25], a meteorological phenomenon that generates intense storms that have resulted in numerous catastrophic floods registered since the XIII century in the Guadalentín basin [2628].

The aquifer recharge comes from the irrigation returns, rainfall, and streams crossing the basin. The basin is located in a low pluviometry area where groundwater inflows are significantly related to seasonal rainfall. Average annual rainfall is 250 mm, with many years with less than 200 mm.

The Guadalentín depression is a fertile basin cultivated since Arabic ages (800 a.C.). In 1960 a wells drilling campaign began in order to improve the area’s productivity by exploiting the important underlying aquifer. Increasing groundwater extractions revitalized Guadalentín agriculture in the 1980s but led to the declaration that the aquifer was temporally overexploited in 1987 [29]. New regulations and water transferred from Tajo’s basin (Tajo-Segura diversion, 1979) produced a reduction of pumping and abandonment of some wells from 1988. Piezometric levels began to slowly stabilize in major areas but continued declining in areas close to the still numerous active wells, concentrated in the North and North-West. From 1960 to 2012 levels declined from near land surface to more than 200 m below land surface. The populations of Lorca and Puerto Lumbreras have continuously increased from 70,000 habitants in 1960 to more than 100,000 in 2012 (data from annual census, Spanish Statistical Office).

Recent studies using SAR from 1992 to 2012 images have revealed a continuous deformation pattern with one of the greatest velocities in Europe (around 11 cm/yr) on its maximum area [15, 17]. The subsiding area is closely related to the thickness of a soft soil layer, presented in the shallower part of the Plio-Quaternary materials and comprising low and very low permeability silts and clays.

As presented, the study area is a populated region with important geological risks and vulnerable structures and infrastructures where a combination of hydrogeological and surface deformation solutions could be used to improve the groundwater management.

3. Data and Methods

3.1. Conceptual Model of Groundwater Flow

During the last decades, the aquifer system declaration as overexploited led water management agencies to compile existing information and generate different versions of the Guadalentín basin numerical groundwater model. The Spanish Geological Survey (IGME, the Spanish acronym) built a numerical groundwater flow model in 1993 covering the period from 1960 to 1992 [30]. The basin management agency, Segura hydrographic confederation (CHS), built its own model in 1991 covering the period from 1982 to 1989, later updated in 2005 with intensive field studies in 2003 and 2004 [31, 32], but with a constant recharge in time. The new model presented in this paper merges all the data from these models, updating and joining them into one model that covers the period from 1960 to 2012.

Following the model build-up, data has been divided into three main groups. The available geological data has been used to define the aquifer system geometry. Recharge and discharge data has been used to analyze the inflows to and outflows from over the aquifer.

3.1.1. Aquifer System Geometry

Geometry of Plio-Quaternary sediments was studied by Cerón [22] through vertical electric soundings (VESs) resulting in the Figure 2(a) thickness map. Recent studies developed by Bonì et al. [15] have demonstrated that Plio-Quaternary materials can be divided into two sublayers. The shallower layer is comprised of soft soils, clays, and silts, while the deeper one is the coarse fraction of Plio-Quaternary materials, sands or gravels, which constitute the main productive layer. Soft soils geometry was estimated by Bonì et al. [15] using geological information from 23 boreholes in the basin (Figure 2(b)). These materials, their location, and differential hydrological behavior have not been used in previous groundwater models developed in the basin. Their inclusion in our model constitutes an important improvement introduced in this study. Although the definition of the basin basement as a horst and graben structure is commonly accepted, the real position is not clearly defined by different authors, with elevations that range from 400 to 1000 m below land surface. In order to simplify the geometry and ensure the stability of the model the lower boundary is located at a medium depth of 500 m below land surface. This third layer is formed by Miocene materials changing from conglomerates and sandstones in the upper zone to marls in the deepest areas.

3.1.2. Inflows

The principal inflows to the aquifer system occur through rainfall surface infiltration and stream infiltration through ramblas. Average annual inflows are summarized in Table 1. According to the rainfall series provided by the Spanish National Meteorological Agency (AEMET) average rainfall is around 250 mm/yr, usually concentrated in storm events from August to October. During the 52 studied years a rainfall declining trend can be observed (Figure 3(b)) with two severe droughts from 1994 to 1996 and from 2001 to 2003. Owing to the presence of an important thickness of clay materials with low permeability, a 30 km2 area in the northern part of the basin has been marked as impervious and not-susceptible to surface recharge (Figure 4(b)).

Another important source of recharge is infiltration from the irrigation return flow. Agricultural importance of this area is linked with water availability, supplied by different sources. This agricultural water demand is being met by water supplied from nearby reservoirs, groundwater extractions, and, more recently, transference of water from other major basins (Tajo-Segura water transference, 1979). Agricultural development has affected the aquifer because of the extractions and the recharge from irrigation returns. Irrigator associations from Lorca (CRL) and Puerto Lumbreras (CRPL) manage the two main irrigated areas of 91.7 and 26.1 km2, located in the northern and southern sectors of the basin, respectively (Figure 4(b)). Irrigation in Puerto Lumbreras increased greatly in the 1960s and 1970s, reaching a peak in 1975. After 1988, nearly related to the restrictions imposed on groundwater extraction, irrigation decreased dramatically, whereas irrigation in Lorca presents a more stable trend during the simulation period. This irrigation area is partially characterized by impervious soils that reduce the useful infiltration surface to 46.8 km2.

Guadalentín seasonal stream and Nogalte rambla, located, respectively, in the north and south of the western margin of the basin (Figure 4(a)), are the most important sources of stream infiltration to the aquifer. The magnitude of these two watersheds has allowed the generation of well-developed hydrogeological connections with the aquifer [14]. There are other minor seasonal streams along the eastern (SE ramblas) and western (SW ramblas) margins that contribute to the aquifer inflows (Figure 4(a)). These streams are too small to be considered as individual entities and to develop important alluvial systems.

3.1.3. Outflows

Discharge from the Alto Guadalentín aquifer occurs through the connection with the Bajo Guadalentín aquifer system (Figure 4(a)) and the numerous wells pumping its resources. Groundwater discharge to streams was not considered because the Guadalentín stream has a losing-disconnected relation with the Alto Guadalentín aquifer system.

Extractions rate (Figure 3(a)) was governed by the increasing irrigation water demand and annual rainfall, reaching its maximum at the end of 1980s when 80 hm3/yr were pumped. After restrictions on groundwater pumpage were imposed in 1989, the extraction rate dropped to less than 30 hm3/yr in the early 1990s. During the mid-to-late 1990s, extractions increased slightly in order to compensate for a series of dry years (Figure 3(b)). From 1960 to 1989 312 wells were registered and distributed throughout the basin (Figure 4(a)). Two areas near the basin center (eastern center and western center, Figure 4(a)) showed denser concentrations of wells. Since 1989 the number of wells has been reduced. Of the remaining 50 wells, the deepest and most productive are mainly concentrated in the western area.

3.2. DInSAR Data

DInSAR processing allows the evaluation of the temporal evolution of surface displacement using SAR images over large areas since the first SAR satellites became operational in 1992 [33]. Considering that piezometric level changes are one of the most important mechanisms that triggered surface deformation [34, 35], the relationship between piezometric level and deformation is being studied due to its potential applications in aquifer resources management [11, 12, 36].

The presented groundwater model recreates piezometric levels history from its original state in 1960 to current 2012 levels. Deformation data do not cover the entire simulated period, but all the periods for which SAR images exist for the study area were processed. Three different sets of images were used from C-band ERS (1992–2000), ENVISAT (2003–2010), and X-band Cosmo-SkyMed (CSK) (2011-2012) satellites.

InSAR-derived ground displacements covering the period 1992–2012 (Figure 5) were collected from previous studies [15, 17, 37]. Deformation measurements from the ERS dataset are from Rigo et al. [16]. This dataset was processed using DORIS interferometric software [38] to complete the coregistration and interferogram generation phases. Time series were computed using the SBAS approach [39] through StaMPS software [40] (Figure 5(a)). CSK dataset (from Bonì et al. [15]) was processed using DIAPASON for the interferometric stage and SPN software during the final products calculation following the Persistent Scatterer approach [41] (Figure 5(c)). Both dataset products, average displacement velocities and line-of-sight (LOS) displacement time series, are described in detail by Gonzalez and Fernández [17] and Bonì et al. [15]. Deformation data from the ENVISAT dataset are from Béjar-Pizarro et al. [37]. This dataset is composed of 27 images from August 2003 to July 2010. It was processed by means of ROI_PAC [42] to focus on the raw data, DORIS to calculate the interferograms and StaMPS SBAS approach to generate interferograms and deformation velocities/time series (Figure 5(b)).

3.3. Groundwater Numerical Model Development and Calibration

In order to develop a groundwater numerical flow model of the Alto Guadalentín aquifer, we have used MODFLOW-2005 with the ModelMuse interface, both developed by the USGS [18, 43]. Modeling was carried out in two steps; first a steady-state simulation was performed using 1960 data, and second, a transient simulation for the period from 1961 to 2012 was performed using the steady-state simulation model as starting conditions. A model time step of one year was specified because the model input data such as water balance, recharge, and extraction were compiled annually.

Model calibration was performed using UCODE_2014 [44] using the graphical interface ModelMate [45] developed by the IGWMC and USGS to minimize the root mean squared residual through adjustment of hydraulic conductivity and storage coefficients.

3.3.1. Model Discretization and Soil Parameters

The model is composed by three layers using square 100 m grid cells over 243.5 km2 across the basin. Flow between cells is controlled by Layer-Property Flow Package (LPF) due to its capability to allow later calibration of hydrogeological parameters. The first layer of the model comprises soft soils Plio-Quaternary materials, the second layer is formed by coarse-grained Plio-Quaternary sediments, and the third layer is formed by Miocene material, as described in Section 3.1.1. Due to the aquifer characteristics and the LPF available options all the layers are described as convertible layers.

Initial hydrogeological parameters have been estimated using data from pumping tests and calibrated information from IGME [30] and CHS [31] models. After the calibration phase is carried out over the presented model, hydraulic conductivity (, cm/s) ranges from to depending on the layer. Coarse-grained Plio-Quaternary materials show higher hydraulic conductivity than Plio-Quaternary soft soils and Miocene layers. The calibrated storage coefficients vary from to with the highest values corresponding to the coarse-grained Plio-Quaternary layer.

3.3.2. Boundary Conditions: Recharge, Discharge, and Initial Conditions

The MODFLOW Recharge Package has been used to approximate rainfall recharge into the aquifer, with an annual average value during the simulated period of 7 hm3/yr and maximum and minimum values from 26.74 hm3/yr and 2.23 hm3/yr, respectively. Rainfall recharge has been uniformly distributed over the 213.5 km2 active upper layer grid cells. Due to the absence of specific infiltration studies we have used data from a previous model implemented by IGME in 1993. In that study a formulation that calculates infiltration as a variable percentage of rainfall between 10% and 20% was adjusted. This percentage is closer to 10% in the dry years and closer to 20% in the rainy years. Using this data from 1960 to 1993, a quadratic regression curve was adjusted and used to calculate the remaining infiltration series from 1994 to 2012.

Irrigation recharge volumes are extremely difficult to acquire and only available in some years. The irrigation volumes for unmeasured years have been interpolated. Infiltrated volumes are calculated as 15% of the irrigated water and introduced in the modeling using the MODFLOW Recharge Package.

A precise estimation of the recharge contribution of Guadalentín and Nogalte stream courses is complex; only one approximation for 1960 [30] exists, being 1.5 hm3/yr. Taking into account the fact that stream flows are generated by the concentration of the rainfall within their basins, it is reasonable to assume that groundwater recharge contributions can be calculated through a direct relation of the pluviometry (Figure 3(b)) and the 1960 estimate. In order to approximate this boundary condition the MODFLOW Flow and Head Boundary Package has been applied [46].

The contribution from western and eastern ramblas is inserted in the model as stream infiltration inflows using the Flow and Head Boundary Package. Western ramblas are located over 20 km between Nogalte rambla and Guadalentín River covering an area of 25 km2. Eastern ramblas reach 38.6 km collecting runoff water from 50 km2. As in the major ramblas, infiltration rate was approximated for 1960, at 1 hm3/yr in the west and 2 hm3/yr in the east [30]. Estimation of the remaining years has been carried out using the rainfall history and a linear relation [30].

The boundary of the Alto Guadalentín and Bajo Guadalentín aquifer system is the natural outflow of the modeled aquifer. Groundwater flows through the shared boundary have been estimated using Darcy’s law and the piezometric data from both sides of the boundary and are simulated in the model using the Flow and Head Boundary Package. Results show that originally water flowed from the Alto Guadalentín aquifer to the Bajo Guadalentín aquifer at a rate of 12 hm3/year. Due to the Alto Guadalentín aquifer overexploitation the original flow decreased and was reversed around 1985.

The main groundwater abstraction is produced by extraction wells disseminated all over the basin. As pumping rates are given as an overall value, estimated rates for the different wells were based on their potential productivity. The northwestern half of the basin presents higher thickness of Plio-Quaternary and Miocene materials with higher permeability than the Southeastern materials, making wells located in the NW deeper and more productive [23, 28]. Using this position factor extraction rates were weighted, allocating higher rates to the northwestern wells. Wells discharge was simulated using the MODFLOW Well Package.

Ideally, a transient groundwater flow model should begin from a steady-state condition. Because pumping from the Alto Guadalentín aquifer began more than 50 years ago, the search for information on steady-state conditions in the aquifer, used as initial condition of the steady-state simulation step, has to go back to the information gathered by the early studies in 1960. That time was chosen because the system was not suffering significant stress. The piezometric level was between 260 and 320 m a.s.l. as seen in Figure 4(c) [29]. Taking into account the elevation of the basin the original piezometric level was from 15 to 60 m depth depending on the area, making it easy to pump with traditional methods.

3.3.3. Calibration Data

Model calibration has been implemented using 47 piezometric level measurement points spatially distributed with 189 annual piezometric measurements recorded (Figure 4(d)) by IGME and CHS during consecutive measurement campaigns. Low measurement/measurement-point ratio and long simulation time indicate that the available time series are not sufficiently long enough to adequately constrain the calibration. Only 4 measuring points have more than 10 years’ data, with the longest one being 23 years. 14 of the remaining points present series from 4 to 10 years of piezometric levels. 29 have piezometric levels with less than 4 years, most of them with only one year. None of the 47 series covers the simulation period. The existing data is well distributed spatially but poorly in time. Additionally near time records are even more limited, bringing to light the necessity of aquifer monitoring improvements. Piezometric level observations were implemented using the Head Observation Package [47].

Piezometric series presented in Figure 4(d) show a continuously declining piezometric level from 1970 to more recent times with a slight recovery in some areas near the late 1980s, consistent with the aquifer history of extractions (Figure 3(a)). Recent measured piezometric levels (2009–2012) revealed a nearly stable trend but with lower levels than older data (average decline of 90 m from 1990 to 2009), indicating a net loss of groundwater storage.

3.4. Deformation Model Development and Validation

Fusion of this groundwater model result with subsidence data resulting from SAR images processing provides the opportunity to analyze both phenomena together and establish relations between them. Using surface deformation, piezometric level history, and soft soil thickness, we propose an empirical formulation, without geotechnical parameters and time independence, which allows the calculation of surface deformation from piezometric level changes. Parameter calibration was carried out using the ENVISAT results because its time centered position (2003–2010) permits the validation of the model fit for both earlier and later time periods using ERS (1992–2000) and CSK (2011-2012) results, respectively. Two independent variables were used to estimate deformation, piezometric level change, and soft soil thickness. Both linear and nonlinear regression models were considered. The quality of the models to replicate the observed data was assessed using coefficient of determination (). Due to the importance of the soft soils layer the percentage of its thickness over complete Plio-Quaternary sediments thickness was calculated and used to distinguish areas with high percentages (>25%) where described analysis and relationships were done. The empirical deformation models were used to estimate deformation within the range of the simulated piezometric level variation. Accumulated deformation in 1992 and 2012 was estimated and validated with ERS and CSK SAR data by evaluating the ratio error/total displacement.

4. Groundwater Model Results and Discussion

Modeling results present the evolution of the piezometric levels throughout 52 years (1960–2012) across the Alto Guadalentín basin and how the actions on the aquifer have affected them. Figure 6 shows the results at specific and relevant dates. In 1972 (Figure 6(a)) during the early expansive phase of the agricultural water use and increased groundwater extractions, piezometric levels were declining slightly and uniformly, and piezometric level depressions began to form in the areas where the wells were concentrated. In 1989 (Figure 6(b)) at the end of the period of high agricultural water use and groundwater extractions when the aquifer was declared partially overexploited, piezometric levels had declined between 100 and 160 m from their original position, creating two steep depressions surrounding the eastern and western well fields. After the measures adopted by the basin management authorities, a decrease in agricultural water use, and consecutive rainy years from 1989 to 1993 (Figure 3), model results in 1993 (Figure 6(c)) reflect a slight recovery in some areas of piezometric levels. The stagnation of the extractive activity in the eastern well field produced the recovery of the piezometric levels over the surrounding area. This fact suggests the ability of the aquifer system to recover its piezometric levels after severe extractions were stopped. Progressive increases in extractions and a lower rainfall period with three severe droughts associated [14] led to declining piezometric levels by the end of the simulation period in 2012 (Figure 6(d)). The active western well field produced a water level decrement of 70 m from 1993 and an overall decrement of 40 m throughout the basin.

In order to evaluate the groundwater model results two important factors must be taken into account: its length, 52 years, and the piezometric level global change, between 170 and 210 m. Both of them can generate important differences between results and observations with slight changes in water contributions and uses (specially, pumping rates and watering). Despite these limitations the computed model levels fit well with the observations. Figure 7 shows the simulated levels against the observed levels at the validation points described in Section 3.1.3. These values group around the line with an RSME of 17.4 m. Average absolute difference between simulated and observed piezometric levels is 14.9 m. Assuming an average piezometric level drop of 190 m, both errors are under 10% of the total water displacement, an acceptable error for a long duration regional model.

An overall comparison of simulated piezometric data reflects a slight overestimation trend in the modeled level (piezometer 2539-1-001) converging both series towards the 1980s (piezometers 2539-2-043, 2539-6-024). From that date the model shows better results as can be seen in piezometer 2539-2-011. Blue and black framed graphics in Figure 8 describe the temporal evolution of the modeled and observed series in some wells, demonstrating model ability to fit well into great piezometric level changes with acceptable errors.

In order to evaluate the spatial distribution of the errors, average errors of each piezometer were interpolated to obtain an error map of the model along the basin (Figure 8) using all the available piezometric data of the 52 modeled years. Highest differences areas are located on the southern edge, a marginal zone without enough data. From center to southwest more important areas also present obvious differences. Piezometric levels simulated at piezometer 2539-1-001 are higher than the observed in the 1973–1980 period (Figure 8). Piezometer 2539-2-041, located in the same area, follows the same pattern during those dates, fitting better after them. Better model fit was obtained in the northern half of the basin reaching errors under 15 or even 10 m along wide areas. This situation can be observed through the 2539-2-011, 2539-2-111, and 2539-2-043 piezometers time series and it is represented in the blue framed graphic (Figure 8).

As seen in Figures 7 and 8, most of the differences are positive, so the model is inclined to overestimate the simulated piezometric levels. This is more noticeable during the first simulation steps (1960–1980) and in the southern half of the basin. Two particular facts could have been underestimated and disturb the model. On the one hand, this area is more arid than the northern half of the basin, its irrigation system is worse developed, and the crops typology adapted to the situation. On the other hand, Puerto Lumbreras irrigated areas grew very quickly during the 1960–1975 period. The combination of these factors can be translated into uncertainty in irrigation returns (less watering for crops adapted to rainfed lands) and extractions (more extractive pressure due to drier conditions and quick development). Both concentrations of wells (eastern and western) are associated with error areas around the model average error (15 m). Piezometers related to the well fields (2539-2-011, 2539-2-043, 2539-3-066, or 2539-3-119) display some of the best fits between modeled and observed series, supporting a precise estimation of water extraction effects in the aquifer system dynamic. The northern area, with denser urbanization and specially affected by subsidence processes [15, 16], is located within the best fitting zone.

5. SAR Deformation: Advance Use of Groundwater Modeling Results and Discussion

Statistical analysis of deformation and piezometric change series results in a low correlation between them (). Despite the poor relationship observed, this study was able to reveal the presence of at least two data populations in the original data.

Previous studies in the Alto Guadalentín basin have already probed the relationship between deformation and soft soil thickness [15, 37]. Taking into account those works and layers geometry described in Figure 2, the percentage of soft soils thickness over the complete Plio-Quaternary materials thickness has been computed. Using this parameter to identify the two observed populations, original SAR deformation data were divided into two samples of 1597 (percentage of soft soils over 25%) and 3393 (percentage under 25%) points each one (Figure 2(c)). SAR deformation points with soft soils percentage over 25% show a stronger relationship () with piezometric changes and were selected to calibrate the deformation model. Note that selected area is one of the lowest error areas in the groundwater model (Figure 8).

Some empirical deformation models were executed during the analysis using water changes as independent variable and, afterwards, introducing soft soils thickness as another independent variable. Using the compressible layer thickness as a variable relies on its strong relationship with SAR-derived deformation, resulting in better fits when using thickness as a linear variable for deformation models, changing from to . In order to improve the contribution of soft soil thickness nonlinear relationships were also considered. Applying thickness as quadratic or exponential variables, determination coefficient () rises to and , respectively. Exponential term for thickness and linear term for piezometric level changes is disclosed as the best formulation to describe the relationship between water, thickness, and deformation (Figure 9). This correlation range of validity is for soft soils thickness over 10 m and piezometric changes from 8.5 m to 25 m. Calculated deformations out of the range are subjected to high error or incoherent results.

An exponential term to describe the water change looks logical watching the data shape, but in this case the formulation would be applied to great piezometric changes (around 190 m between 1960 and 2012), generating deformations extremely high and incoherent. Adopted empirical deformation model is also limited as can be seen in Figure 9(d). Deformations are overestimated with low water changes while high piezometric variations generate underestimated deformations.

Applying the formulations over the ERS validation period (1992–2000) (Figure 10), calculated deformation shows a good spatial relationship with SAR deformation. On the other hand maximum deformation values are underestimated from 28% to 40% in the 8 years’ ERS deformation. Minimum error is associated with the model that explains surface deformation using an addition of a linear correlation with piezometric changes and an exponential correlation with soft soils thickness (Figure 10(c)), supporting it as the most suitable.

Deformation derived from CSK data (2011-2012) was added to the ENVISAT time series due to its short time span and the previously described problems with water changes out of the validity range. Model results for the period 2003–2012 (Figure 11) show differences between the simulated and measured spatial deformations patterns, but in this case simulated maximum values are much closer to measured values (difference under 4%). Spatial discrepancy is mainly related to the influence of the soft soil thickness pattern in the model.

The increasing error from present to past is caused by one of the most important limitations of this model, the temporal position of the calibration period in the deformation history. Material stress due to water withdraw began in 1960 and the calibrated relation calculated during ENVISAT period is related to an inelastic deformation period. Empirical relations calculated are related to inelastic deformations that need less water changes than elastic ones to generate the same deformation and overestimate ancient deformations.

Keeping in mind previously described model limitations and features, deformation accumulated during the period 1960–2012 was calculated. Accumulated deformation in 2012 varies from 3 to 4.5 m. Deformation model was calibrated with a piezometric change validity range where soft soils thickness has a limited influence in the final result also adjusted to that range. If piezometric variations are significantly increased, thickness influence will became low relevant as happened during the period 1960–2003 when piezometric level changed 150 m on average. This effect is much more marked along marginal edges, generating severe deformations in areas with slight movements during monitored periods. Looking for a solution to mitigate this error, whole 52 years’ model results were subdivided into several sets with piezometric level changes into validity range, calculating deformation of each one and joining back the results. Using this alternative way, water changes are in adjustment range improving model results. Enhanced results show an accumulated deformation between 1.1 and 5.4 m during the 1960–2012 period (Figure 12). Although remarkable movements (around 1.5 m) were calculated in near-stable areas, overall deformation pattern is coherent with current observations.

Our results suggest that a maximum deformation of ~5.5 m has occurred since 1960 assuming that model premises and parameters remain constant over time. Marginal deformation is probably overestimated due to deformation model problems in areas to low soft soil thickness. As surface deformation is closely related to soft soil thickness, a better knowledge of the spatial distribution of these materials would improve deformation model results.

6. Conclusions

We have presented here a hydrogeological model of an aquifer system that has experienced intensive pumping activities generating average piezometric level declines of about 150 m during a period of 52 years (1960–2012). The model simulates changes in groundwater flow, from the original steady-state conditions in the aquifer system to current conditions, increasing our knowledge about the aquifer system behavior under stress in order to improve its management. The evolution of piezometry from this model is also a valuable constraint for geomechanical modeling aimed at characterizing the Alto Guadalentín basin that has one of the greatest subsidence rates in Europe (11 cm/yr).

The use of DInSAR has led to a better deformation monitoring due to exploitation of its spatial coverage and high density measurements unattainable by traditional methods. The detailed characterization of the basin materials carried out by Bonì et al. [15] through geological and SAR deformation data analysis has been used to refine the conceptual groundwater model. These geological updates constitute a first contribution of SAR data to improve water management.

Our results also corroborate an important relationship between subsidence, piezometric levels changes and compressible thickness variations. An empirical deformation model based on the groundwater numerical model results was adjusted using the ENVISAT (2003–2010) dataset and then surface deformation during the ERS (1992–2000) and CSK (2011-2012) periods was used to validate this model. In addition, subsidence in the basin over the entire period studied was estimated using the described model ranging from 1 to 5.5 m in 52 years. Estimated subsidence should be understood as maximum values associated with relevant uncertainties derived from groundwater model and deformation model properties, but it constitutes a reference initial point to develop more complex deformational models.

Improvements of the geological data used in this study (e.g., the soft soils thickness map recently updated by Béjar-Pizarro et al. [37]) will allow improvement of our results. Additionally, the coupled groundwater flow and subsidence models, using the Subsidence and Aquifer System Compaction Package (SUB-WT) [48], can be used to refine the original groundwater model and better integrate the deformation data. Thanks to the European Space Agency (ESA) Sentinel-1 constellation, the monitoring of the study area continues today with enhanced characteristics of wide spatial coverage, great temporal resolution (six days’ repeat cycle), and high spatial resolution. Future works combining our results with climatic change scenarios could facilitate the integration of piezometric level and deformation predictions into the management of the aquifer system.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was developed in the framework of AQUARISK Project: Analysis of Geological-Geotechnical Risks due to Groundwater Exploitation Using Space and Terrestrial Techniques. Application on urban structures and infrastructures was funded by Spanish Research Program from Ministry of Economy and Competitiveness (ESP2013-47780-C2-2-R). First author shows gratitude for Ph.D. Student Contract BES-2014-069076.