Quantifying Vertical Deformation in the Tigris–Euphrates Basin Due to the Groundwater Abstraction: Insights from GRACE and Sentinel-1 Satellites

Abstract

This study explores the occurrences of land subsidence in response to dropping groundwater levels in the central part of the Tigris–Euphrates basin. We estimated the groundwater depletion related to human and climate drivers between 2003 and 2017 based on estimates from the Gravity Recovery and Climate Experiment (GRACE) and two global hydrological models (NOAH-3.3 and WGHM-2.2d). The cumulative displacement was calculated using Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) for 96 interferograms, which were generated by 31 images acquired by the Sentinel-1 satellite. The results show that the basin is experiencing a decline in groundwater storage at a rate of −7.56 km³/year with a total loss of 106.81 km³. This depletion rate has led to inelastic compaction and has resulted in subsidence near the city of Baghdad at a rate of −10 mm/year. The measured coherence phase between the two signals is ~0.67, and the depletion precedes the subsidence by ~1.5 months. The new data from GRACE-Follow on, Sentinel-1, and the piezometric water level could help to constrain the rate of depletion and displacements in the basin. Combining these remote sensing techniques provides an independent tool for water management in areas where in-situ data are scarce.

1. Introduction

The Tigris–Euphrates basin (TEB) is transboundary basin bounding two river systems (Tigris and Euphrates), which originate from eastern Turkey, follow through Syria, Iraq, Jordan, Saudi Arabia, and Iran and end in the Persian Gulf [1]. Most of the basin has recently suffered from water risk, especially downstream, which is a result of the extensive use of surface water and groundwater and the frequent occurrence of droughts (Figure 1). For example, the dry spell from 2007 to 2009 caused declines in flow and reservoir storage to less than 70% of their normal levels [2,3].

The scarcity of surface water is putting more pressure on groundwater aquifers, which has resulted in a higher extraction rate [4] and decreases in storage by −25.6 ± 3.7 km³/year, with a total loss of −76.9 ± 10.1 km³ [3]. In theory, the decline in groundwater level to historical levels without replenishment could lead the aquifer sediments to compact and the land to subside [5]. Such phenomena are persistent in urban and agriculture regions throughout the world, such as the San Joaquin Valley [6], Santa Clara Valley, California [7], Mekong Delta [8], México [9], and northern Beijing Plains [10].

In the surrounding regions of the basin, such as Iran, the subsidence related to the drying of groundwater wells has been a persistent problem in the last decade [11]. In the Rafsanjan plains, the subsidence rate is estimated to be >5 cm/year and has increased to 30 cm/year in local regions based on the interferometric synthetic aperture radar (InSAR) technique [11]. This research investigates whether the aquifers in the TEB are suffering inelastic compaction due to the massive groundwater extraction, the decline of water levels, and the progressive drought in recent years.

We address this question by first providing evidence for the groundwater depletion and land subsidence in the TEB. Then, we prove that such subsidence is related to the rapid extraction of the groundwater in recent years. The impetus for this investigation is the scarcity of real in-situ data, including piezometric water levels and sufficient subsidence monitoring systems (e.g., levelling systems and the Global Positioning System (GPS)). To overcome this issue, we combine information from the geodetic satellites (Gravity Recovery and Climate Experiment (GRACE) and Sentinel-1) with Global Hydrological Models (GHMs) to better understand the changes in groundwater storage (GWS) and explore the climate and human-driven depletion in the basin.

GRACE is sensitive to the total mass variation based on temporal changes of the gravity field [12]. Thus, it can provide the changes and trends of GWS with the help of GHMs. Sentinel-1 satellites offer a tool to estimate the crustal deformations at the millimeter scale, which can be seen as a change of the hydraulic head in aquifer regions. We used the Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) [13] for data from the Sentinel-1 satellite between 2015 and 2018 to resolve the temporal and spatial subsidies in the central part of TEB. To determine the relationship between the depletion and land subsidence, we filtered the two-time series in time and space, and we found the common power, coherence, and consistent phase by applying the cross-wavelet transform [14].

2. Materials and Methods

2.1. Estimating the Groundwater Storage (GWS) Change

GRACE satellites have been used over the last 16 years to investigate the changes of GWS indirectly with the help of GHMs, such as the Global Land Assimilation System (GLDAS). GRACE represents the integrated vertical column of the water change, including the soil moisture ($\Delta \mathrm{SMS}$), snow water equivalent ($\Delta \mathrm{SWE}$), surface water storage $\left(\Delta \mathrm{SWS}\right)$, canopy ($\mathrm{CW}$), and GWS, while GLDAS simulates these storages except for GWS and SWS. GWS can be obtained by subtracting the GLDAS simulated storages from the total water storage (TWS) obtained with GRACE based on the water budget approach [16].

$$\Delta {\mathrm{TWS}}_{\overline{\mathrm{GHM}}}=\left(\Delta \mathrm{SMS}+\Delta \mathrm{SWE}+\Delta \mathrm{CW}+\Delta \mathrm{SWS}\right)$$

(1)

$$\Delta \mathrm{GWS}=\Delta {\mathrm{TWS}}_{\overline{\mathrm{GRACE}}}-\Delta {\mathrm{TWS}}_{\overline{\mathrm{GHM}}}$$

(2)

The errors of the GWS based on Equations (1) and (2), can be formulated as

$$\sigma \mathrm{GWS}=\sqrt{\begin{array}{c}{\left(\sigma \Delta {\mathrm{TWS}}_{\overline{\mathrm{GRACE}}}\right)}^2+{\left(\sigma \Delta {\mathrm{TWS}}_{\overline{\mathrm{GHM}}}\right)}^2\end{array}}$$

(3)

where $\sigma$ is one sigma errors of the storage.

We used the GRACE mascon solutions, which use the mass concentrations as a prior constraint to derive the reduction of the temporal changes of the gravity field and increase the signal-to-noise ratio [17,18,19]. To represent the TWS in the TEB between January 2003 and June 2017, we used the average of the three mascon solutions released by three different institutes: the Center for Space Research, the University of Texas at Austin (UT-CSR), NASA Jet Propulsion Laboratory (JPL), and NASA’s Goddard Space Flight Center (GSFC). For other components, we used the outputs of the GLDAS-2.1 (NOAH-3.3) model for the same period. The NOAH-3.3 model provides optimized states for the storages at a depth of 2 m for four soil layers [20] but does not simulate $\Delta \mathrm{SWS}$ (e.g., rivers, lakes, wetlands, and reservoirs). Therefore, we obtained other outputs from the WaterGAP Global Hydrology Model (WGHM 2.2d). WGHM uses water balance estimation (e.g., irrigation, water management) to simulate all storages and water fluxes, such as surface and subsurface runoffs, river discharge, lakes, wetlands, GWS, human withdrawals, snow, rivers, and canopy) [21,22,23]. The data are available at a resolution of 0.5° (55 km at the equator) for the period of 1979 to 2016. We used two versions of the model and considered human impact to emphasize the role of human interventions in the basin (e.g., irrigation and groundwater extraction).

2.2. Temporal Land Deformation in Tigris–Euphrates Basin (TEB)

The response of the aquifer system to the decline of the water level was investigated using the coherence small baseline subset [24] and atmospheric corrections for the Sentinel-1 data. Coherence SBAS relies on the combination of redundant series of unwrapped interferograms with good coherence that can be attained using pairs with narrow temporal and spatial baselines. The deformation signal $\Delta \tau$ in the interferometric phase $\Delta \varphi$ at time $i, j$ is contaminated by phase errors $\alpha$ (e.g., atmospheric and ionospheric errors, antenna noise), and other errors $ϵ$ arise from sources such as the topography according to Equation (4) [25,26].

$$\Delta {\varphi }_{ij}=\Delta {\tau }_{ij}+{\alpha }_j-{\alpha }_i+{ϵ}_{ij}$$

(4)

Assuming that $\alpha$ has a random distribution, the errors can be estimated by applying common point stacking [26] and then subtracting from the interferometric phase, which leads to increases in accuracy and strengthens the deformation phase. We used 31 image datasets from the C-band of the Sentinel-1 satellites, which are operated by the European Space Agency (ESA) [27]. The datasets represent an ascending orbit with 72 paths and 106 frames acquired from 26 May 2015 to 10 April 2018 over the central part of the TEB (Figure 1 and

Table 1. List of Sentinel-1A images used in this study. PBL is the parallel baseline, and PPBL is the perpendicular baseline of the satellite in meter.

Date	PBL (m)	PPBL (m)
2015.05.26	0.00	0.00
2015.06.07	19.41	6.72
2015.06.19	71.91	105.48
2015.07.13	72.54	124.01
2015.07.25	15.69	−27.02
2015.08.18	45.56	75.30
2015.11.22	71.82	109.82
2015.12.16	66.46	107.44
2016.01.09	76.53	160.95
2016.02.02	51.28	102.53
2016.03.21	58.83	75.62
2016.04.14	45.73	10.76
2016.05.20	125.45	145.82
2016.06.13	41.98	32.23
2016.07.07	78.92	70.49
2016.09.17	120.16	131.20
2016.09.29	93.99	111.05
2016.10.23	10.86	31.45
2016.11.04	30.30	71.45
2017.11.29	36.82	96.68
2017.12.11	19.12	77.74
2017.12.23	15.12	69.71
2018.01.04	49.60	125.45
2018.01.16	67.80	143.60
2018.01.28	74.42	151.14
2018.02.09	44.62	103.60
2018.02.21	45.66	97.30
2018.03.05	31.60	76.26
2018.03.17	75.01	131.62
2018.03.29	88.12	145.53
2018.04.10	38.63	52.36

). The SAR images are co-registered pixel by pixel based on a geometric approach while accounting for topographic variations. This is done with the help of a 30 m digital elevation model (DEM) [28] and the provided orbital information [25,29].

The 31 SAR images are aligned to a super master image captured on 11 April 2016, which has the minimum difference in baseline elevations between the images (Figure 2). The 96 interference patterns between the images are constructed based on the minimum spatial baseline of 100-m in the acquired data. The topographic phase and the flat-earth effect are removed by using the 30-m resolution DEM. Next, a Gaussian low-pass filter is applied with 300 wavelengths. The unwrapped phase is obtained by using the Statistical-Cost Network-Flow Algorithm for Phase Unwrapping (SNAPHU) software [30]. We generated a rate deformation map that represents the cumulative vertical displacements in the TEB from 2015 to 2018 based on stacking the unwrapped phases and the time series of the deformations at the acquisitions dates. The time and frequency phase of the mean series of vertical displacements are compared with the derived GWS depletion using a cross-wavelet transform to identify the possible connection with ground deformations in the TEB.

2.3. Cross-Wavelet Transform between Depletion and Displacements

The continuous-wavelet transform (CWT) is a band-pass filter technique to extract the possible physical relationships between two geophysical time series in time-frequency space based on the findings of a consistent phase angle (time lag), common power, and coherence [14]. CWT can be defined as:

$$D\left(\frac{\left|W_n^X\left(s\right)W_n^{Y*}\left(s\right)\right|}{\sigma X\sigma Y}<p\right)=\frac{Z_v\left(p\right)}{v}\sqrt{P_K^X{P}_K^X}$$

(5)

where $W_n^X$ and $W_n^Y$ are continuous wavelets that represent the relative phase between the two time series in X and Y (the derived GWS and the vertical displacements) with length $n$ and standard deviation $\sigma .$ $Z_v$ is the level of confidence associated with the probability (p) density function of the two product root square distributions, and P_K is the background power spectrum.

The correlation is significant between the two signals when the left-hand side of Equation (5) is larger than the maximum power spectrum of the red noise. The cross-wavelet transform (XWT) treats the Morlet wavelet as a mother function, which balances between localizing the signal in time and frequency. The extracted significant periodic signal is represented by the cone of influence (COI), where the phase discontinuity drops by $e^{-2}$ at the boundary of the COI. The CWT is used to define the relationship between different geophysical processes—for example, the Baltic maximum sea ice extent and Arctic Oscillation index [14], land subsidence and GWS volume change [31], and climate, vegetation factors, and ΔTWS [32]. We employed the CWT to expand the groundwater volume changes and vertical displacements with corresponding time into a time-frequency space. Then, the XWT is employed based on the two CWTs to expose the regions in time-frequency space with a local correlation locked phase under a statistical significance level. The statistical significance was determined by the Monte Carlo method [14].

3. Results and Discussion

Figure 3 shows the time series of the storage components and the derived GWS. GWS is plotted with the uncertainty obtained from Equation (3). We also present the simulated GWS from WGHM-2.2d under two scenarios regarding the impact of human intervention. The long-term trend was calculated using the Sen slope [33] under a 95% confidence interval and p-value < 0.0001. From January 2003 to June 2017, the long-term rate of the GWS depletion in the TEB basin is about −9.93 ± 2 mm/year, which is equal to −7.56 km³/year with a total loss of 106.81 km³. Much of this loss occurred after 2007, when the basin experienced an intense drought.

The mean lower uncertainty for the GRACE-GWS is −39.71 mm, and the upper uncertainty is approximately −22.30 mm. The spatial variations of the TWS and GWS trends show that the maximum depletion occurs in central section of the basin as opposed to the northern and southern sections (Figure 4). The depletion of TWS has a range of −11 to −21 m/year in the central section and decreases to <−7 mm/year in the upper stream in Turkey and to <−9 mm/year downstream in Iraq. Similarly, the GWS has a depletion rate of 11 to −21 m/year in the central and downstream sections of the basin. The upper stream has a relatively small rate of <−7 mm/year.

The WGHM estimate of GWS was compared the result derived from the GRACE, which indicated that the model underestimates the component magnitude. However, it still indicates the long-term trend and the change of the GWS compared with GRACE-GWS. The measured correlation between the GRACE-GWS and WGHM simulated GWS is R = 0.84 for GWS-HI (with human interventions) and R = 0.55 for GWS-NHI (without human interventions). One can infer that groundwater depletion is occurring in the TEB and is exacerbated by the drought conditions and the rapid human extraction that occurred in the last decade (2007–2017).

The displacement rate obtained by applying the SBAS to the SAR images ranges from −40 to 40 mm/year. Higher subsidence was detected near Baghdad with a rate of about > –10 mm/year, as well as the eastern part of the study area (Figure 5). The maximum subsidence is located in superficial urban regions, wetlands, and paddy fields, where the overexploitation of groundwater is occurring. The lithology of the subsided area is represented by quaternary sediments (e.g., sand, silt, and clay), and the groundwater depth ranges from 2 to 50 m [34]. The flow of channels near this region has decreased markedly since the 1960s and dropped to 520 m³/s by 2012 (Supplementary Materials, Figure S1) [35,36]. For this reason, the extraction of GWS and the dependence of irrigation systems have increased.

Such withdrawal of groundwater could be mainly based on groundwater during drought periods when surface water resources decrease. For example, following 2007–2009, the Iraqi government and citizens dug more than 1000 new wells. In 2018–2019, when the drought became worse, it is more likely that such activities on groundwater have increased. In general, the withdrawal of groundwater in recent years has resulted in Iraq ranking as 21th among countries with the highest groundwater depletion in the world in terms of withdrawal ratio to the water supply, which is >80% [37]. The withdrawal of surface water and groundwater constitutes more than 97.2% of the total withdrawal in Iraq [4].

One possible interpretation of the vertical displacements in TEB is the oil and gas production, which can induce such land displacements. In fact, frame 1 covers the giant oil field of East Baghdad. However, no vertical displacement has been reported in the area because of oil production. In light of the absence of information about the production rate of the field, there is a possibility of the contribution of oil extraction to the land subsidence near Baghdad. However, the same anomaly propagates into frame 2, where there is no oil or gas production. Therefore, we can infer that the anomaly is most likely caused by the rapid extraction of groundwater and the compacting of aquifer materials. A pronounced uplift of about 28–30 mm/year was observed in the northern area of frame 1 (non-urban region) and the eastern area of frame 2 (the western sections of Zagros Mountains).

The results of the continuous and cross-wavelet analyses are presented in Figure 6. As expected, the dominant significant period is the seasonal component (12 months) with a power spectrum of 6–8 over the time series of 2003–2017 for the GWS changes. For the land displacement, the significant period is 1.5 to 2 months from November 2015 to July 2016 with the same power spectrum range of 6–8. The common power between the two signals is in a period of 1.5 to 2-months. The coherence of the two signals in frequency space is about 0.67 for the 1.5–2 months, which is above the 5% confidence level.

The direction arrows are in the down-and-left direction, which implies that the two signals are out of phase and that the first signal (GWS changes) leads the second (displacements). The phase angle is about ~1.5 months. The results from GRACE-TWS show that the trend of the TWS from 2013 to 2016 has a decreasing rate of −21.09 mm/year, while the GWS has a negative rate of −2.41 mm/year from January 2011 to December 2016. In the same period, subsidence occurs, and the temporal patterns of the two signals are coherent. However, due to the short time of the comparison (i.e., 2015–2017), such comparisons could be misleading.

Based on the spatiotemporal results of groundwater depletion in TEB, the subsidence in urban regions, and the coherence level of the depletion and land displacements, we conclude that the subsidence in the center of the TEB has a mutual relationship with the activities of groundwater extraction.

However, some limitations are worth mentioning. First, we compared the GWS depletion across the whole basin with local subsidence, which can be seen as a rough assessment due to the coarse resolution of GRACE (>200,000 km²) and the short term of comparison. Future work must include the results of GRACE-FO with finer resolutions and adopt many Sentinel-1 frames to convert the displacements to hydrological information.

Second, obtaining a continuous deformation rate is challenging and is mainly based on the availability of the data. In this analysis, the deformation rate resulted from two different interferogram networks, and the resulting rate is not unique. The quality of the interferogram networks involves a tradeoff between the temporal and spatial baseline information with respect to the temporal changes of phenomena (e.g., groundwater seasonal changes). Lastly, even though the wavelet analysis does not rely heavily on the length of the data since the periodicity is the dominant component, it is safer to use a long time series of data.

4. Conclusions and Implications

We examined the subsidence induced by groundwater depletion in the central TEB using remote sensing observations and GHMs. The TEB experienced a GWS decline rate of −7.56 km³/year with a total loss of 106.81 km³ from 2003 to 2017. The GWS depletion induced by human intervention surpasses the impact of climate on the TEB. The resulting depletion has induced a subsidence rate >−10 mm/year near Baghdad. The depletion and subsidence signals covary with coherence values of approximately 0.67, and the depletion precedes the subsidence by less than two months. The new data from the GRACE-FO as well as updated Sentinel time series and piezometric data from in-situ water levels could help to constrain the rate of depletion and displacements in basins. These remote sensing techniques using GRACE and Sentinel still provide independent tools to practice water management in places where in-situ data are scarce.