Understanding the dynamic behaviour for the Madrid aquifer (Spain): insights from the integration of A-DInSAR and 3-D groundwater flow and geomechanical models

Advanced Differential Interferometric Synthetic Aperture Radar (A-DInSAR) techniques and 3-D groundwater flow and geomechanical models are integrated to improve our knowledge about the Tertiary detritic aquifer of Madrid (TDAM). In particular, the attention is focused on the Manzanares-Jarama well field, located to the northwest of Madrid, which experienced five cycles of extensive groundwater withdrawal followed by natural recovery, to cope with the droughts occurred in summer 1995, 1999, 2002, 2006, and 2009. Piezometric records and A-DInSAR data acquired by ERS-1/2 and ENVISAT satellites during the periods 1992–2000 and 2002–2010, respectively, have been used to calibrate the groundwater flow and the geomechanical models. A time-lag of about one month between the hydraulic head changes and the displacements of the land surface has been detected by a joint wavelet analysis of A-DInSAR and piezometer head time series. Overall, the results show the effectiveness of the proposed integrated approach composed of A-DInSAR and 3-D geomechanical model to characterize the aquifer-system response during and after the groundwater withdrawal.

1 Introduction

Groundwater numerical models provide a quantitative framework to understand the hydrogeological processes and represent a practical tool to comprehend how an aquifer system works (Anderson et al., 2015). Groundwater models integrate numerically the Partial differential equations governing the physical process of the water movement within a defined portion of a porous medium. Modelling design requires to set-up the model boundaries and geometry, provide the aquifer parameters and choose the code able to simulate the physical processes of interest (Anderson et al., 2015). Therefore, the groundwater model conceptualization needs extensive datasets that often lack due to the high costs required to collect.

Figure 1 Geographical location from Google Earth, Data SIO, NOAA, U.S., Navy, NGA, GEBCO, US Dept of State Geographer, ©2019 Google, Image Landsat/Copernicus (a) and geological map (b) of the study area. The piezometers and the well locations are also depicted. The geological map is derived from IGME (1995) and the water bodies and the rivers from https://servicio.mapama.gob.es/sia/visualizacion/descargas/mapas.jsp (last access: 2 March 2020).

To remediate, at least in part, this paucity of measurements, in this work a novel methodology using A-DInSAR data as supporting tool for groundwater modelling purposes is presented. The methodology has been developed and tested to understand the dynamic behaviour of the Madrid aquifer. It represents a strategic groundwater resource in Spain, exploited in the past to cope with an increased groundwater demand during drought periods. The achieved results give insight about a fast response of the aquifer-system during and after the extraction periods.

2 Study area

The study area is the Manzanares-Jarama wellfield, that is located northwest of Madrid, in Spain (Fig. 1). The analysed wellfield includes 32 wells, which withdraw groundwater from the “Manzanares-Jarama aquifer”, and 17 observation wells. The analysed aquifer extends for about 540 km² and represents only a portion of the 2500 km² Tertiary detritic aquifer (TDAM) located in the tectonic depression of the Madrid basin. The aquifer is a heterogeneous multi-layer system characterized by sand lenses embedded in a clay and clay-sand matrix (Ezquerro et al., 2014). The aquifer includes Tertiary detritical deposits and a Quaternary coverage (Martínez-Santos et al., 2010). The northern and southern boundaries of the aquifer are the Paleozoic relief of the Sierra of Guadarrama and evaporitic materials, respectively (Hernández-García and Custodio, 2004). Conversely, in the western and eastern boundaries areas, the alluvial deposits of the Manzanares and Jarama rivers are found, respectively.

Figure 2 Line-of-sight (LOS) velocity maps measured usingERS-1/2 (1992–2000) (a) and ENVISAT (2003–2010) images (b). World imagery Source: Esri, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

3 Data and methodology

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 3 Piezometric head time series acquired in the period 1998–2010. The location of the observation wells is reported in Fig. 1b.

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 4 Schematic overview of the integrated approach used to reconstruct the 3-D model with the support of A-DInSAR data.

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):

$$\begin{array}{}\text{(1)}& S_k={\displaystyle \frac{\mathrm{\Delta }d}{\mathrm{\Delta }h}},\end{array}$$

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):

$$\begin{array}{}\text{(2)}& \mathrm{\nabla }\left(k\cdot \mathrm{\nabla }h\right)=S_{\mathrm{s}}\cdot {\displaystyle \frac{\partial h}{\partial t}}+q\end{array}$$

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):

$$\begin{array}{}\text{(3)}& G{\mathrm{\nabla }}^{\mathrm{2}}{u}_i+\left(\mathit{\lambda }+G\right){\displaystyle \frac{\partial \mathit{\epsilon }}{\partial i}}={\displaystyle \frac{\partial p}{\partial i}}i=x,y\end{array}$$

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:

$$\begin{array}{}\text{(4)}& G={\displaystyle \frac{E}{\mathrm{2}\left(\mathrm{1}+\mathit{\nu }\right)}}\end{array}$$

and

$$\begin{array}{}\text{(5)}& E={\displaystyle \frac{\left(\mathrm{1}+\mathit{\nu }\right)\left(\mathrm{1}-\mathrm{2}\mathit{\nu }\right)}{c_{\mathrm{b}}\left(\mathrm{1}-\mathit{\nu }\right)}}\end{array}$$

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 Results and discussion

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 5 3-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.

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.

5 Conclusions

In this work, the potential of the integrated use of A-DInSAR data and 3-D groundwater flow and geomechanical models to capture and assess aquifer dynamics is shown. The approach has been applied to investigate the response during and after pumping of a portion of the Madrid aquifer. The aquifer-system has been characterized by a time-lag of about one month only between the hydraulic head changes and the land displacements, suggesting that the clayey units within the aquifer-system play a minor role.