3.1 Satellite and piezometric data analysis

The analysed satellite data are Synthetic Aperture Radar (SAR) images acquired by ERS-1/2 and ENVISAT in descending and ascending mode, respectively. The ERS-1/2 images cover the period from April 1992 to November 2000, and the ENVISAT dataset between August 2003 and September 2010. The ERS-1/2 and ENVISAT scenes were processed using the PSP-IFSAR technique (Costantini et al., 2008) to obtain the A-DInSAR dataset including line-of-sight (LOS) displacement time series and average velocity for each target in the periods 1992–2000 and 2003–2010 (Fig. 2).

Available observation well measurements recorded the hydraulic head changes from 1998 to 2012 (Fig. 3) and show three major extraction phases causing a groundwater drawdown of about 100–150 m, corresponding with major drought periods followed by recovery periods (Yélamos and Villarroya, 1991; Canal de Isabel II Gestión, 2014).

Figure 3Piezometric head time series acquired in the period 1998–2010. The location of the observation wells is reported in Fig. 1b.

Download

The methodology implemented by Tomás et al. (2016), based on a wavelet approach, was used for the joint analysis of the A-DInSAR and piezometric level time series to understand cause-and-effect relationships and to identify intermittent periodicity. On the other hand, Continuous Wavelet Transform (CWT), was applied to identify nonstationary intermittent oscillations in time series data. Cross Wavelet Transform (XWT) and Wavelet Coherence (WTC) were applied to find the common power and relative phase in time-frequency space of two-time series (Grinsted et al., 2004). The result is a 2-D plot of complex numbers. A high absolute value suggests a common time patterns between both time series. Furthermore, the phase (i.e. the angle of the vectors) of the XWT provides a quantification of the time shift (delay time) between the two time-series. It is worth nothing that both the A-DInSAR and piezometric head time series have been preliminary equalized in the time sampling interval using a 35 d interval that corresponds to the revisit time of ERS and ENVISAT. Then, a linear regression of the time series records was used to distinguish linear and non-linear component. The wavelet approach has been applied using the non-linear component resulting from the difference between the original time-series and the fitted linear component.

Figure 4Schematic overview of the integrated approach used to reconstruct the 3-D model with the support of A-DInSAR data.

Download

3.2 Integration of the A-DInSAR into numerical models

In this work, the A-DInSAR data are used as supporting tool for different phases of the numerical model reconstruction (Fig. 4). In particular, the A-DInSAR data are exploited to (i) define the lateral boundaries of the static model, (ii) estimate the model parametrization, (iii) choose an appropriate code to simulate the physical processes, and (iv) to calibrate of the model. The joint use of the wellfield location and the average velocities measured by the A-DInSAR technique allows to identifying the influence area of the wellfield. Furthermore, the piezometric head time series and the displacement time series have been jointly used to estimate the aquifer storage (S_s) using a methodology applied by several authors on A-DInSAR-based and extensometers-based displacement measurements (Sneed and Galloway, 2000; Burbey, 2003; Tomás et al., 2010; Bonì et al., 2016). More precisely, the specific storage is computed dividing the bulk storage (S_k) by the aquifer thickness. The bulk storage is estimated using the following equation Eq. (1):

where Δd is the displacement, estimated by the satellite data, and related to the aquifer compaction as a result of the hydraulic head decrease Δh. The joint analysis of A-DInSAR displacements and piezometric head time series, as described in Sect. 3.1, supports the choice of the simulator in the geomechanical model.

Finally, A-DInSAR displacement time series are used for cross-comparison with the simulated displacement to verify the model results.

3.3 Numerical model approaches

The numerical approaches implemented in this work aims simulating groundwater flow and pumping induced land displacements. First of all, a static model is build up to reconstruct the three dimensional geometry including the main hydrogeological layers of the study area. Then, the 3-D static model is used for the hydrogeological and the geomechanical simulations. More precisely, the 3-D finite element (FE) grid of the selected domain is constructed using GEN3-D (Teatini et al., 2006). The groundwater flow model is simulated through the SAT3-D (SATurated Three-Dimensional) code (Teatini et al., 2011) that solves the following equation Eq. (2) (Bear 1972):

where k and h are the hydraulic permeability tensor and the hydraulic head, respectively, t is time, and q is the source/sink term. S_s is the specific storage, which is related to oedometric bulk compressibility c_b through the relationship $S_{\mathrm{s}}=\mathit{\gamma }\left(c_{\mathrm{b}}+\mathit{\varphi }\mathit{\beta }\right)$ with ϕ the medium porosity and β the volumetric water compressibility. The output of the groundwater flow model is the incremental pore pressure variation (p) due to groundwater pumping that is used as input in the geomechanical model. By assuming an elastic isotropic medium, the SUB3-D code (Teatini et al., 2006) is used to estimate the incremental displacement field u along the x, y, and z directions. SUB3-D solves the following equations Eq. (3) (Verruijt, 1969):

where λ and G are the Lamé constant and shear modulus of porous medium, respectively, $\mathit{\epsilon }=\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}$ is the volumetric strain, and u_i is the displacement component along the ith coordinate direction.

It is worth noting that:

and

where E and ν stand for the Young modulus and the Poisson ratio, respectively.

The groundwater flow and geomechanical models are calibrated using a trial and error approach. The hydraulic conductivity was calibrated the records at the piezometric wells starting from values reported by Hernández Torrego (2012) and from values obtained via pumping test for pilot points (e.g. CB5 piezometer, see the location in Fig. 1). The Young modulus and Lamè constants were calibrated using InSAR outcome. Consequently, the solution of Eq. (1) must takes the elastic storage obtained using Eqs. (4) and (5) starting from literature values (Ezquerro et al., 2014). Furthermore, for the geomechanical model the Poisson's ratio ν was set to 0.3, representing a common value used in the case of sedimentary deposits (Teatini et al., 2006).

4.1 Remote assessment of aquifer-system response during and after pumping

The analysis of the A-DInSAR displacement and the piezometric head time series has been performed for 17 piezometers of the Manzanares-Jarama wellfield using the wavelet approach as described in Sect. 3.1. The results provide insight about a mean time lag of 1.0±0.6 months, with minimum and maximum values equal to 0.0 and 3.3 months, respectively. The spatial pattern of the aquifer-system response is slightly larger in the northern sector of the study area. The aquifer-system response during and after pumping is fast, showing land subsidence during the extraction phase coinciding with the drought periods, and uplift in the recovery phase, providing fundamental information about the governing parameters controlling the hydro-geomechanical response.

Figure 53-D model of the Manzanares-Jarama wellfield (a) and simulated and observed displacement at piezometer FE-1. (b) See the location in Fig. 1. The observed displacements are an average value of the values measured on the PS located within a buffer zone of 1000 m from the well, therefore also the standard deviation is reported.

Download

4.2 Numerical simulation of aquifer-system response during and after pumping

The join analysis of the satellite and piezometric data has provided fundamental information to design and calibrate the 3-D model of the Manzanares-Jarama well field. In particular, the model has been built including three main hydrogeological units; i.e. the upper phreatic aquifer, the intermediate aquitard and the lower confined aquifer (Fig. 5a). The model domain extends for the upper 700 m of the sedimentary sequence, whose top and bottom maps have been extracted by the ground surface elevation using the high resolution 5 m cell digital elevation model from Instituto Geográfico Nacional (IGN) and geological information (IGME, 1989). Since the piezometric drawdown caused by the wellfield did not reach a distance larger than 10 km from the wells, (Ezquerro et al., 2014; Béjar-Pizarro et al., 2017), the horizontally domain extends in order to account a buffer ring of 10 km around the wellfield. The aquifer's dynamic has been simulated during the period from 1994 to 2010, that includes four 1–1.5 year long periods of groundwater withdrawal. The modelling results show that the periods of groundwater extraction caused a land subsidence higher than 20 mm over an area of about 190 km² and a maximum horizontal movement up to 6 mm. The simulated displacement generally agrees with the A-DInSAR measurements, capturing the minimum and maximum peak values (Fig. 5b). Furthermore, the calibrated model has allowed quantifying the land surface movements in the central zone of the study areas where A-DInSAR displacements are unavailable due to the presence of vegetation coverage.