Interpolation of GPS and Geological Data Using InSAR Deformation Maps: Method and Application to Land Subsidence in the Alto Guadalentín Aquifer (SE Spain)

Abstract

Land subsidence resulting from groundwater extractions is a global phenomenon adversely affecting many regions worldwide. Understanding the governing processes and mitigating associated hazards require knowing the spatial distribution of the implicated factors (piezometric levels, lithology, ground deformation), usually only known at discrete locations. Here, we propose a methodology based on the Kriging with External Drift (KED) approach to interpolate sparse point measurements of variables influencing land subsidence using high density InSAR measurements. In our study, located in the Alto Guadalentín basin, SE Spain, these variables are GPS vertical velocities and the thickness of compressible soils. First, we estimate InSAR and GPS rates of subsidence covering the periods 2003–2010 and 2004–2013, respectively. Then, we apply the KED method to the discrete variables. The resulting continuous GPS velocity map shows maximum subsidence rates of 13 cm/year in the center of the basin, in agreement with previous studies. The compressible deposits thickness map is significantly improved. We also test the coherence of Sentinel-1 data in the study region and evaluate the applicability of this methodology with the new satellite, which will improve the monitoring of aquifer-related subsidence and the mapping of variables governing this phenomenon.

1. Introduction

Groundwater is a strategic resource that represents the primary source of fresh water for many communities worldwide [1]. Groundwater is stored in aquifer systems that can suffer compaction from groundwater pumping, especially unconsolidated alluvial aquifer systems, causing land subsidence [2]. This subsidence can be accompanied by numerous damaging effects, for instance fluvial and coastal flooding, earth fissuring, damage to buildings and infrastructures and groundwater storage loss [3,4,5,6]. Identifying and mapping the different variables involved in the subsidence phenomena is crucial to understand, model and make predictions about the evolution of the ground deformation.

Groundwater numerical models can be used to reproduce the behavior of aquifers, but in order to obtain realistic results, a dense mapping of the variables involved in the phenomena is required [7]. In situ measurements, such as GPS displacements, piezometric data and geologic boreholes, are very valuable, but expensive and can be difficult to obtain. Consequently, when they exist over a region, in situ measurements are usually very sparse.

Over the past two decades, Interferometric Synthetic Aperture Radar (InSAR) techniques have revolutionized the study of ground displacements associated with different phenomena, including aquifer system compaction and ground subsidence [8,9]. InSAR has been used to detect and monitor many areas of subsidence produced by groundwater exploitation [10,11,12], deduce aquifer hydraulic properties [13,14,15] and estimate hydraulic head evolution [16].

Geostatistical techniques are an extension of linear regression models to estimate georeferenced variables at unsampled locations. The common geostatistical interpolation approach, known as kriging, takes into account the spatial structure through variogram models at the same time, which provides unbiased estimates with minimum variance. One extension is Kriging with External Drift (KED), which uses a linear regression model with another (external) georeferenced variable for the trend. The KED method has been applied to rainfall estimation [17], piezometric contouring [18], temperature mapping [19] or hydraulic conductivity estimation [20] using different variables as the external drift, including elevation, groundwater modelling solutions or electrical resistivity. However, this method has never been applied using InSAR data as the external drift.

In this paper, we apply the KED method to interpolate GPS vertical deformation rates and the thickness of compressible deposits in the Alto Guadalentín aquifer (SE Spain). This region is affected by the fastest subsidence associated with groundwater pumping in Europe, with deformation rates higher than 10 cm/year. We first describe the study area geological location and the main characteristics of the Alto Guadalentín aquifer (Section 2). After summarizing the data and methods used (Section 3), we apply the KED approach to our two datasets and obtain continuous maps of the two variables (Section 4). Finally, we compare the new maps with previous studies and discuss the future application of the method in light of the excellent Sentinel-1 coherence and wide spatial coverage (Section 5).

2. Study Area

The Guadalentín basin, SE Spain, is one of the driest regions of Europe [21]. This basin is a tectonic depression located in the eastern part of the Betic Cordillera, an alpine orogenic belt resulting from the collision between the African and Iberian plates. The convergence of ~5 mm/year [22] between the two plates is accommodated by large NE-SW faults, such as those controlling the mountain fronts limiting the basin [23,24]. The tectonic rate of the large faults in the study region (e.g., The Alhama de Murcia fault, [24]) is ≤1 mm/year.

The substratum of the basin consists of Permian-Triassic metamorphic materials exhibiting a horst and graben structure [25]. The metamorphic substratum is covered by materials of Miocene age that include marls, sands, conglomerates and calcarenites. On the top of these, there is a succession formed by conglomerates, sands, silts and clay constituting Plio-Quaternary alluvial fans [26]. The Alto Guadalentín detritic aquifer, located between the cities of Lorca and Puerto Lumbreras, developed in the Plio-Quaternary deposit. The aquifer system has a surface of 275 km², and its thickness varies between 100 and 300 m ([27]; Figure 1). Based on stratigraphic information from 23 boreholes, [28] classified the Plio-Quaternary materials filling the basin into compressible and non-compressible, the former consisting of clay and silt layers and the latter formed by conglomerates and sands (see Figure 1c). The continuous pumping of groundwater in the aquifer, mainly for agricultural use, led to a decrease in the piezometric levels of more than 200 m since 1975 [28], resulting in the overexploitation of the aquifer and the land subsidence studied here.

Several studies have focused on deciphering the factors controlling this phenomenon by analyzing land subsidence [28,29,30]. The most comprehensive study to date [28] analyzed 20 years of deformation time series (1992–2012) in combination with groundwater level data from 1975–2012 and geological information. These authors did not find a clear correlation between groundwater drawdown and ground subsidence during the period 1992–2012. Instead, they found a good correlation between groundwater drawdown during the previous period (1975–1992) and cumulative displacements from 1992–2012, so they concluded that the displacements observed from 1992–2012 was triggered by the piezometric drawdown that began in the 1970s, suggesting a long delayed compaction process. Their comparison of the magnitude of the displacement rate and the thickness of the Plio-Quaternary sediments filling the basin revealed that there is not a direct relationship between them, which they attributed to the spatial variability of the grain size and the thickness of the Plio-Quaternary deposits. Finally, Boni et al. [28] compared the thickness of the compressible soils (within the Plio-Quaternary deposit) with the magnitude of the deformation rate and found a clear correlation between them. They concluded that the aquifer system is experiencing a slow consolidation process that was triggered by the piezometric drawdown that started in the 1970s and is controlled by the presence of a very thick layer formed by low permeability and compressible deposits. This behavior of the aquifer system, with groundwater level variations triggering the subsidence, whereas local geological conditions control the magnitude of the deformation, has been observed elsewhere (e.g., [31]).

Our study is based on the well-established relationship between the thickness of the compressible deposit and the ground deformation found by [28]. We first confirm the existence of this correlation with our particular InSAR dataset, and then, we use a different interpolation method to improve the compressible deposit thickness map created by [28].

3. Data and Methods

3.1. Estimation of InSAR Deformation Rates

Ground deformation in the study area was measured using ENVISAT SAR data from [32]. The ENVISAT SAR sensor (C-band) has an average incidence angle of 23°; the spatial resolution is 25 m; and its minimum revisit time is 35 days [33]. Twenty seven ENVISAT ASAR images acquired from ascending Track 373 between 2003 and 2010 were processed using a small baseline approach [34]. We did not process the descending track covering the study region because acquisitions between 2006 and 2010 are very scarce. We first focused the raw data using ROI_PAC [35] and then, 74 interferograms were produced using DORIS [36]. The topographic phase contribution was corrected in the interferograms using the Shuttle Radar Topography Mission data [37]. The selection of stable scatterer pixels and the time series analysis was performed using StaMPS [38]. From this processing we extracted the average Line-Of-Sight (LOS) subsidence rate maps covering the two areas used in the interpolation: the area labelled A in Figure 1, covered by 35,875 points, and the area labelled B in Figure 1, covered by 24,584 points. We used two overlapping areas, only covering the area with the sparse point measurements to be interpolated in each case, to optimize the calculation time during each interpolation.

3.2. Estimation of GPS Deformation Rates

The study area is crossed by four lines of the high precision levelling network of Spain (Red Española de Nivelación de Alta Precisión or REDNAP) managed by the Spanish National Geographic Institute (Instituto Geográfico Nacional de España or IGNE). The orthometric height of 39 points of this network located in the study region were measured by the IGNE during three different precise levelling campaigns: 20 points were measured in August–September 2004 (red diamonds in Figure 1b); 11 points were measured in March 2005 (yellow diamonds in Figure 1b); and 8 points were measured in March 2009 (blue diamonds in Figure 1b). The exact date of each measurement is indicated in Supplementary Table S1. The resulting heights have mm precision (1.5 mm × $\surd K$, where ‘K’ is the levelling line length expressed in kilometers).

We measured these 39 levelled benchmarks using the GPS rapid static surveying method between 28 January and 1 March 2013. All of the measurements were made using a Leica System 1200 GNSS receiver and an ATX1230 antenna with a 1-Hz logging rate. All sites were occupied for a minimum period of 10 min. This allows obtaining an uncertainty in the vertical relative position of 2–3 cm. In addition, two continuous GNSS stations (LORC and LOR1; Figure 1) were used as reference sites. These stations were selected by their proximity to the study area (less than 20 km), which provides a maximum uncertainty of 2 cm in the ellipsoidal height calculation. All of the ellipsoidal heights are referred to the Geodetic Reference System 1980 (GRS80) ellipsoid.

The reference vertical positions provided by the IGNE and the new vertical positions measured in 2013 were compared to estimate the vertical displacement rate of each site. For this purpose, we converted the ellipsoidal heights into orthometric heights using the official EGM2008-REDNAP geoid model [39] with an absolute uncertainty of 3.8 cm. Considering the maximum uncertainty associated with the ellipsoidal heights, we obtain a maximum uncertainty of 4.3 cm in our orthometric heights, which is sufficient considering the high displacement rates measured in the study region [28].

3.3. Compressible Deposit Thickness

The thickness of compressible deposits was compiled in [28] using stratigraphic information from 23 boreholes. For each borehole, the total thickness of silt and clay layers located above the layer of non-compressible materials, mainly formed by gravels, was estimated.

For this study, 5 out of the 23 boreholes were removed from the analysis: two of them because of their location outside the Alto Guadalentín aquifer and the other three because they presented low confidence value (due to low quality description of the borehole). The 18 boreholes selected for this study are shown in Figure 1. They present a good spatial distribution, covering the complete area of the Alto Guadalentín aquifer. Among the 18 boreholes, 8 of them have 1 m of compressible deposit thickness, and 2 of them have 0 m. The other 8 boreholes have compressible deposit thickness values ranging between 4 and 192 m (Supplementary Figure S1). The distribution of the compressible deposit across the basin and its relationship with the non-compressible deposit and the substratum is shown in Figure 1c.

3.4. Kriging with External Drift

To create continuous maps from our discrete variables (GPS vertical velocities and compressible deposit thicknesses), we applied the kriging with external drift (KED) method. Kriging is a method of interpolation that minimizes the error variance using a weighted linear combination of the data [40]. KED performs the standard kriging algorithm considering that the drift is locally represented by some auxiliary variable (InSAR deformation rates in our study). This variable needs to be known everywhere in the field, or at least to be sampled densely enough so that the value at any point can be obtained consistently. This is our case, in which an ordinary kriging has been applied to the 35,875 InSAR data to obtain a complete InSAR image.

As in any other kriging procedure, a model of the spatial correlation is required. In the KED, this model has to be inferred knowing that the InSAR information serves as a local drift. The model has been calculated in the scope of the intrinsic random functions of order k theory [41], using polynomial isotropic generalized covariances. A nugget effect was avoided to honor the GPS information during the interpolation. To apply KED, the fundamental hypothesis of a linear relation between the two variables must make physical sense.

The geostatistical study was performed using the Isatis geostatistics software [42]. For each of the discrete variables (GPS vertical velocities and compressible deposit thickness), the analysis followed these steps: First the linear relation between InSAR data and the discrete variable is verified. Then InSAR data are interpolated to obtain a continuous map of deformation velocity, using the ordinary kriging method. Finally, the kriging with external drift is applied to two discrete variables (GPS deformation rates and compressible deposit thickness) using a continuous map of InSAR-derived deformation as the external drift.

To validate the interpolated maps, the “leave-one-out” or cross-validation method is applied [43]. It estimates data “blindly”, removing one datum at a time, and predicting the associated data where values are known. The predicted (Z*) and real (Z) values at the location of the omitted point are compared. The kriging standard deviation estimation (S*) represents the error predicted by the model and depends on the model and the location of the neighboring information [44]. Dividing the cross-validation error by S*, the standardized error is obtained, allowing one to compare the magnitudes of both the real and the predicted error [45]. All of these quantities can be plotted to show how good the prediction is. If the prediction errors are unbiased, the mean prediction standardized error should be near zero. Then, if the standardized error lies outside a given interval, the point requires some attention. Here, we define the interval [−2.5; 2.5] to focus on the 1% extreme values of a normal distribution. Another key value to evaluate the quality of the interpolation is the variance of the standardized error, which should be close to 1. This means that the experimental error is close to the model error, and thus, the model parameters are suitable [46].

An additional tool to analyze the error data is to test for the “conditional unbiasedness”. To check for conditional bias, we plot the errors as a function of the interpolated values. Ideally, this scatter should result in no correlation and no increase in the estimated variance [47].

4. Results

4.1. InSAR and GPS Deformation Rates

Figure 2 shows the InSAR-derived map of deformation rates, measured along the satellite line-of-sight, GPS vertical deformation rates and profiles for the Alto Guadalentín basin. The InSAR rate map shows a deformation pattern elongated in the 45°N direction, compatible with the deformation field measured in previous studies [28,29,30]. Maximum deformation rates (11 cm/year away from the satellite) are located in the center of the Alto Guadalentín basin (black star in Figure 2a and LOS time series in Figure 2e). Near the areas with maximum deformation rates, there is a lower density of points due to decorrelation caused by very large displacements. A different processing approach (e.g., the joint exploitation of phase and amplitude [48]) could have improved the spatial density of measurements in those regions.

InSAR measures changes in the line-of-sight distance between the ground and the satellite (23° from the vertical for ENVISAT). To compare these measurements with GPS vertical rates, we transform LOS into vertical displacement. Considering that groundwater subsidence mostly produces vertical deformation and that InSAR data are more sensitive to vertical than horizontal deformation, we assume that the main contribution to the displacement field is from vertical deformation. This is supported by the InSAR and GPS deformation rates measured over the two permanent GPS stations LORC and LOR1 (Figure 1a) by [28]. Both datasets cover the same temporal intervals (2007–2010 and 2011–2012), and their results showed that deformation rates were very similar, with average difference <5 mm/year. We convert LOS into vertical displacement by dividing the InSAR measurements by the cosine of the ENVISAT incidence angle (23°).

Vertical deformation rates estimated from GPS data are represented by blue arrows in Figure 2a. As explained in Section 3.2, these deformation rates are obtained by comparing the reference vertical positions provided by the IGNE, which varies in the network between August−September 2004 and March 2009 (see Figure 1b), and the new vertical positions measured in 2013. Values vary between 0 cm/year and −12 cm/year (Supplementary Table S1), and their distribution is very similar to the InSAR deformation field (Figure 2).

Figure 2b–d shows the comparison of vertical deformation rates from GPS and InSAR data. In general, vertical rates from GPS data are higher than vertical rates from InSAR. Note that InSAR and GPS data do not cover the same temporal period, which could explain this difference. Our vertical velocities, obtained using the GPS rapid static mode, confirm the subsidence detected with InSAR. The use of this method is acceptable in this case due to the high entity of the deformation phenomenon, but it is not applicable elsewhere. Specifically, when calibrating InSAR data, the use of data from permanent GPS stations is more appropriate.

Note that, as indicated in Section 2, the tectonic rate of large faults in the study region is very low (≤1 mm/year) and consists mainly of horizontal displacement. Therefore, although the tectonic deformation is contained in our data, it is undetectable.

4.2. Correlation InSAR-GPS and InSAR-Compressible Deposit Thickness

Figure 3a shows the correlation between InSAR and GPS velocities at the 39 GPS points measured. Figure 3b shows the comparison between InSAR and compressible soil thickness at the 18 points where the thickness of compressible materials was compiled from boreholes.

InSAR and GPS velocities present a high correlation, with a Pearson coefficient value of 0.96, which indicates that the fundamental hypothesis of a linear relation between the two variables is met, and therefore, the KED can be applied. The correlation of InSAR and GPS measurements makes physical sense because both measure ground deformation. Even though the component of the deformation field and the time interval is different in each case, the high correlation suggests that ground deformation is dominated by vertical subsidence, and its spatial distribution does not change significantly over the two periods.

A high correlation between InSAR and compressible deposit thickness is also apparent in Figure 3b when considering the whole dataset (Pearson coefficient value of −0.94), but it is less clear for points with very low values of compressible deposit thickness. Specifically, among the 18 boreholes selected for the interpolation, eight of them have 1 m of soft soil thickness, and two of them have 0 m (see Supplementary Figure S1). InSAR-derived deformation over these 10 points with low compressible deposit thickness varies between −0.2 and −3.1 cm/year. If we do not consider the point with velocity −3.1 cm/year (G-115 in Supplementary Figure S1), the rest of the points are located in regions with deformation rates between 0 and 1 cm/year, which is within the stable range (±1 cm/year) proposed by [28]. Therefore, for nine of the 10 points with low compressible deposit thickness, deformation rates can be considered stable.

The correlation of InSAR-derived deformation and the thickness of compressible alluvial deposits was well established by [28] after analyzing groundwater level data, InSAR-derived deformation and detailed geological information of the basin. Our analysis confirms the existence of this correlation with our InSAR dataset.

4.3. Map of GPS Vertical Deformation Rates

The continuous GPS velocity map interpolated from 39 GPS data using the KED method is shown in Figure 4. Maximum subsidence rates (13 cm/year) are located in the center of the Alto Guadalentín basin. The interpolated values have been compared with vertical deformation rates from InSAR and original GPS vertical deformation rates in the cross-section (Figure 4b–d).

As previously noted when we compared InSAR and GPS vertical rates in Section 4.1, GPS deformation rates are generally higher than InSAR deformation rates. For instance, in Cross-section A-B (Figure 4b), a difference between InSAR and GPS deformation rates is observed for a distance greater than 17 km: while InSAR data in this region indicate no significant ground deformation, GPS measurements show subsidence rates of ~1 cm/year. This different deformation rate could be explained by a deformation phenomenon affecting the region during the period 2011–2013, only covered by the GPS data. However, since the GPS deformation rates are within the error bar of GPS measurements, the difference between both deformation rates could be also attributable to the GPS uncertainty.

The cross-validation results for the GPS interpolation are shown in Figure 5. The distribution of standardized errors is represented in Figure 5a. The errors range between −1.85 and 2.24, with no values outside the interval −2.5–2.5, indicating that the interpolation procedure correctly re-estimates the data values. The mean error is close to zero (−0.02 cm/year), and the variance of the standardized error is close to one (1.03), indicating the good quality of the interpolation

Figure 5b shows the errors as a function of the interpolated GPS rates to check for the conditional unbiasedness. No correlation exists between the errors and the GPS rates, and no increase in estimation variance is observed. Therefore, no conditional bias occurs in our interpolation.

4.4. Compressible Deposit Thickness Map

Figure 6a shows the resulting map of compressible deposit thickness in the Alto Guadalentín aquifer. Maximum thickness values (greater than 100 m) are concentrated in the center of the Alto Guadalentín basin, where the maximum deformation rates are located. This map can be compared with the map in Figure 6b, created by [28] using a very similar dataset, but a different interpolation approach that did not take into account the distribution of InSAR measurements. The difference between both maps is shown in Supplementary Figure S2. These results are discussed in Section 5.2.

Figure 7 shows the results of the cross-validation procedure for the interpolation of the compressible deposit thickness. The histogram in Figure 7a shows the distribution of standardized errors, which does not appear to follow the normal distribution. This can be attributed to the significant number of points with close-to-zero values of compressible deposit thickness (Supplementary Figure S1). Nevertheless, errors range between −2.01 and 1.92, with no values outside the interval −2.5–2.5 indicating that the interpolation procedure correctly re-estimates the data values. The mean error is close to zero (−0.053 cm/year), and the variance of the standardized error is close to one (1.16), indicating the good quality of the interpolation.

Figure 7b shows the distribution of interpolated compressible deposit thickness versus the standardized error to check for the conditional unbiasedness. Different from the interpolation of GPS vertical rates, here, we observe some correlation for compressible thickness values lower than 50 m. This conditional bias can be attributed to the small number of boreholes available in the interpolation. Although in this case, the quality of the interpolation is not as good as the interpolation of the GPS rates, we obtain a global unbiased map, with a mean error close to zero, thanks to the KED method and the InSAR data.

5. Discussion

5.1. Comparison of Vertical Deformation Rates with Previous Studies

In the last few years, several studies have focused on the Alto Guadalentín basin subsidence. Here, we compare our map of vertical ground deformation rates with previous deformation measurements.

González and Fernández [29] studied deformation rates in the Guadalentín basin using ERS and ENVISAT SAR measurements from 1992–2007. Their InSAR deformation maps showed a very similar pattern to our deformation maps. They detected maximum deformation rates during the period 1996–2001 (~15 cm/year) and a decreasing rate during the period 2003–2007, with maximum values of ~10 cm/year, in agreement with our InSAR and GPS values.

Rigó et al. [30] computed nine differential interferograms during the period 2004–2005 using ENVISAT data. They obtained the same pattern of deformation as this study, but slower deformation rates (maximum value −7 cm/year). This difference can be explained by the different time interval (2004–2005 versus 2003–2010) and the different processing approach (differential interferograms versus time series) used in [29].

More recently, Boni et al. [28] studied the deformation in the Alto Guadalentín basin using InSAR and GPS data. They estimated vertical velocities of two permanent GPS stations in the region, LOR1 and LORC (Figure 1). For the LORC station, the velocity estimated for the period 26 August 2007–22 May 2013 is −9 cm/year [28], which is slower than the velocity estimated in our GPS map for the same location (−12 cm/year). For the LOR1 station, the vertical velocity estimated for the period 4 September 2008–15 June 2013 is −1 cm/year [28], which is also slower than the vertical velocity in our map for the same point (−3 cm/year). These differences can be due to the different time interval in [28] and this study (2004–2013). InSAR measurements from [28] cover the period 1992–2012 and include SAR data from ERS and ENVISAT satellites for the period 1992–2007 (from [29]), the ALOS satellite for the period 2007–2010 and the COSMO-SkyMed (CSK) satellite for the period 2011–2012. In Figure 8, we compare our InSAR and GPS results with InSAR deformation rates for periods 2007–2010 and 2011–2012 from [28]. The ALOS satellite (L-band) has an average incidence angle of 34°, 10-m spatial resolution and a minimum revisit time of 46 days [49]. The CSK satellite (X-band) has an average incidence angle of 42°, ~5-m spatial resolution and a minimum revisit time of one day [50]. LOS displacements are converted into vertical displacement by dividing the InSAR measurements by the cosine of the incidence angle of each satellite. GPS deformation rates for the period 2004–2013 are in close agreement with ENVISAT deformation rates covering the period 2003–2010. However, the ALOS data show a decreased deformation rate during the period 2007–2010 that continues in the period 2011–2012 covered by COSMO-SkyMed data. This deceleration in deformation rates cannot be detected with our GPS dataset due to its low temporal resolution. It was however detected by the permanent GPS station LORC, which is located near the region of maximum deformation (Figure 1a), which measured a change in deformation rates from 10.4 cm/year in the period 2007–2011 to 7.9 cm/year in the period 2011–2012 [28]. This deceleration in deformation rates has been suggested to be caused by a delayed consolidation process triggered by the intense groundwater drawdown (215 m) that occurred in the aquifer between 1975 and 1992 [28].

5.2. Comparison of the Compressible Deposit Thickness Map with Previous Studies

The cause of the high subsidence rates detected in the Alto Guadalentín aquifer seems to be a combination of two factors: the groundwater drawdown in the aquifer between 1975 and 1992, which triggered the phenomenon, and the distribution of the compressible deposit, which controls the magnitude of the deformation [28]. The location of the fine-grained sediments, thus, can be used to predict where the higher rates of deformation will occur. Conversely, in this study, we use the detected deformation rates to infer the distribution of the compressible deposit thickness over the basin.

Figure 6b shows the compressible deposit thickness map estimated by [28]. Some of the differences observed between the two maps can be attributed to the different dataset used in both studies: while the map in Figure 6b was built using the 23 available boreholes for the interpolation, in this study, we removed data from five boreholes because of their low quality and their location outside the Alto Guadalentín aquifer. However, most of the differences between the two maps are due to the interpolation approach used in each case: while we used the KED approach, [28] used an ordinary kriging interpolation method to create the map in Figure 6b and, therefore, did not take into account the InSAR information. Using the KED approach, the interpolation takes into account not only the information at 18 points, but also the spatial distribution of more than 20,000 points in the InSAR deformation map. Supplementary Figure S2 shows the difference between both maps. Positive residuals have an average value of 20.5 m and are located along the longitudinal axis of the Alto Guadalentín basin, while negative values, with an average value of 12.6 m, are concentrated in the northeastern part of the aquifer.

5.3. Future Research Directions: Sentinel-1 Data

The applicability of the methodology presented here is subjected to two conditions: (1) the availability of InSAR data covering the study region with high coherence; (2) the availability of measurements of some discrete variable highly correlated with InSAR data. In this study, we use GPS and compressible soil thickness data as the discrete variable, because they are highly correlated with InSAR in the study region. Piezometric data are also available in the Alto Guadalentín aquifer, but the relationship with ground deformation is not clear [28].

The availability of InSAR data with high coherence depends on the study region and the characteristics of the available SAR data (band, revisit time). SAR data from European C-band ERS-1/2 and ENVISAT satellites provide global coverage of the land surface for the period 1992–2010, but due to the long revisit time of these satellites (35 days), temporal decorrelation can occur, resulting in low coherency within the deformation maps [51].

The Sentinel-1 satellites, successors of the ERS-1/2 and ENVISAT satellites, with a 12-day repeat orbit cycle (six days when using both Sentinel-1A and Sentinel-1B satellites) and small orbital baselines provide greater spatial coverage and temporal sampling than the previous European missions [52,53,54]. The shorter revisit time improves the coherence over the same area with respect to other sensors [51,55]. Our preliminary result from Sentinel-1A data (Supplementary Figure S3) demonstrates the excellent coherence achieved over a wide area with the new satellite, a key condition to apply the proposed methodology.

6. Conclusions

Subsidence due to groundwater exploitation affects many aquifers in the world, and future climatic conditions combined with population growth are expected to aggravate this situation. Improvement in the characterization of the variables implicated in the phenomena requires both more measurements and better strategies to interpolate them.

Exploiting the high measurement density and the spatial coverage of InSAR deformation maps, this study illustrates how the KED method can be applied to map the spatial distribution of highly correlated discrete measurements. In our cases study, these variables are GPS vertical velocities (39 measurement points) and compressible soil thickness values obtained from 18 boreholes over the Alto Guadalentín aquifer in SE Spain.

Our first result provides a continuous GPS velocity map for the Alto Guadalentín basin, which is in agreement with previous measurements in the region, indicating that vertical deformation during the period 2004–2013 was concentrated in the center of the basin. Even if the low temporal resolution of our GPS measurements does not allow detecting the gradual change in deformation rates identified with InSAR and GPS data by [28], this result shows how the combination of discrete GPS and InSAR measurements can improve monitoring strategies in areas affected by wide subsidence. Our compressible soil thickness map improves the previous map from [28] and constitutes an excellent input for groundwater numerical models.

The preliminary InSAR results from Sentinel-1A data suggest that the methodology proposed here will be applied more easily and robustly with the new satellite, which will allow one to obtain better coherence over wider regions thanks to its great spatial coverage (250 km width in the Interferometric Wide Swath or IW mode), fast repeat cycle (six days) and high spatial resolution (5 m × 20 m for the IW mode). In this context, the combination of discrete sparse measurements of variables related to groundwater extraction with InSAR deformation maps will improve the monitoring of land subsidence due to aquifer exploitation and mapping of discrete variables governing this hazard.