Satellite Monitoring of Mass Changes and Ground Subsidence in Sudan’s Oil Fields Using GRACE and Sentinel-1 Data

: Monitoring environmental hazards, owing to natural and anthropogenic causes, is an important issue, which requires proper data, models, and cross-validation of the results. The geodetic satellite missions, for example, the Gravity Recovery and Climate Experiment (GRACE) and Sentinel-1, are very useful in this respect. GRACE missions are dedicated to modeling the temporal variations of the Earth’s gravity field and mass transportation in the Earth’s surface, whereas Sentinel-1 collects synthetic aperture radar (SAR) data, which enables us to measure the ground movements accurately. Extraction of large volumes of water and oil decreases the reservoir pressure and form compaction and, consequently, land subsidence occurs, which can be analyzed by both GRACE and Sentinel-1 data. In this paper, large-scale groundwater storage (GWS) changes are studied using the GRACE monthly gravity field models together with different hydrological models over the major oil reservoirs in Sudan, that is, Heglig, Bamboo, Neem, Diffra, and Unity-area oil fields. Then, we correlate the results with the available oil wells production data for the period of 2003–2012. In addition, using the only freely available Sentinel-1 data, collected between November 2015 and April 2019, the ground surface deformation associated with this oil and water depletion is studied. Owing to the lack of terrestrial geodetic monitoring data in Sudan, the use of GRACE and Sentinel-1 satellite data is very valuable to monitor water and oil storage changes and their associated land subsidence over our region of interest. Our results show that there is a significant correlation between the GRACE-based GWS anomalies (ΔGWS) and extracted oil and water volumes. The trend of ΔGWS changes due to water and oil depletion ranged from –18.5 ± 6.3 to –6.2 ± 1.3 mm/year using the CSR GRACE monthly solutions and the best tested hydrological model in this study. Moreover, our Sentinel-1 SAR data analysis using the persistent scatterer interferometry (PSI) method shows a high rate of subsidence, that is, –24.5 ± 0.85, –23.8 ± 0.96, –14.2 ± 0.85, and –6 ± 0.88 mm/year over Heglig, Neem, Diffra, and Unity-area oil fields, respectively. The results of this study can help us to control the integrity and safety of operations and infrastructure in that region, as well as to study the groundwater/oil storage behavior.


Introduction
Assessment of surface mass changes is valuable for natural resources management, hazard preparedness, and food security. Long-term monitoring of water and oil depletions using spacebased data is important for natural resources and hazard management. Gravity Recovery and Climate Experiment (GRACE) data have great potential to estimate total water storage (TWS) and are used comprehensively to study hydrological processes [1,2] and groundwater storage (GWS) changes due to, for example, water depletion [2,3], permafrost thawing [4,5], and vertical land motion [6]. Mass transportation can contribute to vertical land motion, that is, uplift or subsidence. Land subsidence is generally caused by a various set of natural and human causes, for example, mining, extraction of liquids near the Earth surface (e.g., petroleum extraction), permafrost thawing, and seismic activities. However, in some cases, natural processes generate relatively slow vertical movements, for example, land uplift in Fennoscandia and Laurentia [6], where, on the contrary, subsidence due to human activities is relatively rapid [7]. The classical techniques of subsidence monitoring, for example, repeated precise leveling, are extremely time-consuming and expensive, especially in large areas, in addition to their low spatial resolution. However, in contrast, the interferometric synthetic aperture radar (InSAR) technique has been widely used to measure the ground surface deformation caused by different geological and geophysical processes [8]. Development of the time-series algorithms and its positive impact of reducing atmospheric artifacts [9,10] enhanced the contribution of the InSAR to land subsidence associated with groundwater depletion.
A key characteristic affecting most environmental studies in Africa is the lack of observations, particularly about sustained time-series data, which gives a considerable rise to uncertainty in the regional statistics. Excessive mass extraction (e.g., water) may result in problems related to water resources and environmental issues such as water resources conflicts and, in some cases, ground subsidence [11,12]. Numerous studies have been carried out recently on groundwater and hydrological aspects using different remote sensing techniques and in situ data. For example, Castellazzi et al. [13] investigated the GRACE/InSAR combination procedure for Central Mexico, where huge groundwater depletion was reported. They concluded that the most notable uncertainty lies in the selection of the GRACE data processing approach. They also presented InSAR-derived aquifer compaction maps. Rahaman et al. [14] proposed, developed, and tested a model to forecast the GRACE-derived groundwater anomalies for sustainable groundwater planning and management in the Colorado River, where a reasonable numerical results are achieved. Bonsor et al. [15] focused on the Nile basin and the use of the groundwater recharge model to interpret the seasonal variations in the terrestrial water storage indicated by the GRACE. Fallatah et al. [16] showed a substantial groundwater depletion that was attributed to both climatic and anthropogenic factors by quantifying the temporal variation of freshwater resources in the arid Arabian Peninsula in addition to identifying factors affecting these resources. Rateb and Kuo [17] explored the occurrences of land subsidence in response to groundwater level changes in the central part of the Tigris-Euphrates basin using GRACE and SAR data. Depletion was computed at a rate of 7.56 km 3 /yr for the GWS, which resulted in a subsidence near the city of Baghdad at a rate of 10 mm/yr. Becker et al. [18] jointly analysed the GRACE, satellite altimetry data, and precipitation, in order to study the spatiotemporal variability of hydrological parameters over the East African Great Lakes region. They found that the TWS change from GRACE and precipitation displays a common mode of variability at an interannual time scale, with a minimum in late 2005, followed by a rise in 2006-2007, which was attributed to the strong 2006 Indian Ocean Dipole (IOD) on East African rainfall. They also showed that the GRACEbased TWS was linked to the El Niño-Southern cycle and the combination of the altimetry data with the TWS allows estimating soil moisture and groundwater volume variations over the lakes drainage basins. Bonsor et al. [19] used GRACE monthly solution provided by three different processing centres to examine the TWS changes in 12 African sedimentary aquifers, to study relationships between the TWS and rainfall, and to estimate the GWS changes using four hydrological models, for example, the variable infiltration capacity model (VIC). They found that there were no continuous long-term decreasing trends in the GWS from 2002 to 2016 in any of the African basins. Du et al. [20] reported their findings based on SAR and GRACE data for monitoring ground surface subsidence owing to groundwater extraction and underground mining activities in the Ordos Basin, China. Ouma et al. [21] integrated GRACE measurements and Global Land Data Assimilation System (GLDAS) models in order to study the inter-annual variations and the GWS changes in the Nzoia River Basin in Kenya. They found that GRACE and GLDAS models could provide reliable data sets suitable for the study of small to large basin GWS variations, especially in areas with limited in situ data. Ahmed et al. [22] analyzed GRACE solutions over North and Central Africa along with relevant data sets (e.g., topography) to assess the utility of GRACE data for monitoring elements of hydrologic systems on local scales. They found that the hydrologic cycle is largely controlled by the estimated mass changes and GRACE data can be used to investigate the temporal local responses of a much larger scale of hydrologic systems.
Land subsidence due to oil extraction has also been studied using SAR data. Xu et al. [23] observed 15 cm subsidence in the center of the Lost Hills field Southern California during a 105-day period. They concluded that this surface subsidence is because of the vertical shrinkage of the reservoir, which resulted in the pore pressure drop. Fielding et al. [24] reported subsidence rates as high as 40 mm in 35 days or >400 mm/yr owing to oil extraction in the San Joaquin Valley, California. Tamburini et al. [25] investigated surface deformation by the persistent scatterer InSAR (PSInSAR) method for monitoring of surface deformation in reservoir area because of fluid extraction from, or injection into, subsurface reservoirs. In most cases, oil field subsidence may not cause a large impact on surface structures as underground mining performs. However, the long-term effect of oil fields subsidence on the local ecological environment, the subsurface hydrological system, and the surface drainage pattern can be significant [26].
The objective of this study is to use available observations from different satellite remote sensing techniques to study the solid Earth response to mass transportation owing to current climate change and oil/water extraction in Sudan. Hereinafter, the term Sudan refers to the territory that includes Sudan and South Sudan. The satellite-based data are collected from low-orbited dedicated missions, which provide a homogeneous global coverage. We study the rate and extent of mass change and its corresponding ground subsidence using GRACE and Sentinel-1 SAR data. For this purpose, the mass variations are estimated in terms of the GWS change and surface deformation in the study region. Firstly, large-scale TWS anomalies (ΔTWS) are studied using GRACE monthly solutions [27,28]. However, GRACE data suffer from inherent data processing problems and need to be addressed properly to achieve a reliable and accurate product. As satellite gravimetry is not sensitive to the Earth's geocenter motion, degree one coefficients, which represent the position of the Earth's instantaneous center of mass relative to an Earth-fixed reference frame, need to be added to the GRACE products by a model of the degree one variations [29]. Independent estimation of geocenter motion that is not quantified in the GRACE solutions is a source of uncertainty in GRACE-based mass change estimation. Moreover, comparing the time series of fluctuations in C20 coefficients observed by GRACE and those achieved by an analysis of satellite laser ranging (SLR) tracking data shows that there are significant differences over time. Therefore, as SLR-based C20 coefficients are more accurate, GRACE-based C20 coefficients should be replaced by those introduced by SLR solutions [30]. Accordingly, the uncertainty of the second degree of zonal spherical harmonic coefficients is another possible source of uncertainty in GRACE-based mass change estimation. Furthermore, as GRACEderived Stokes coefficients are contaminated by spatially correlated errors [31], especially at smaller wavelengths in which the errors are significantly larger than the signal, surface mass change estimates from GRACE monthly solutions are corrupted by north-south stripe noise [32]. Hence, to filter out the effect of the correlated noises on the surface mass change estimates, one needs to utilize the Stokes coefficients that are decorrelated by applying appropriate filters. In this study, we used different isotropic and non-isotropic filters to filter out the effect of the correlated noises in the monthly-normalized spherical harmonic coefficients. We utilized both non-isotropic DDK filters [33,34] and a conventional Gaussian spatial smoothing function [35,27] with 300 and 500 km smoothing radiuses to investigate the effect of applying different filters on our GWS estimates over regional scales. Another task to be tackled is studying different hydrological models such as GLDAS and the WaterGAP Global Hydrology Model (WGHM) [36,37] for the estimation of ΔGWS.
Secondly, the permanent scatterer interferometry (PSI) technique [38] is used to estimate the surface ground deformation rate owing to water/oil depletion and drought in some selected oil production reservoirs in Sudan (see Figure 1). The long record of oil production data, collected from Jan 2003 to Sep 2012, is also utilized in this study.
The paper is organized as follows. The study area is briefly summarized in Section 1. The data, such as GRACE, oil wells production record, and SAR data, and the processing methods are presented in Section 2. The results and discussions are provided in Section 3. On the basis of the results, the relation between the oil extraction and the estimated GRACE-based ΔGWS and then the ground subsidence is shown. This will be valuable information as one can use it for improving future hydrological models over the areas with heavy oil productions. Furthermore, these types of studies enrich our knowledge of reservoir behavior and help achieve more efficient reservoir management and prediction of future performance with obvious economic benefits.

Study Area
The country of Sudan is approximately 29% desert, 19 % semi-desert, 27% low rainfall savanna, 14% high rainfall savanna, 10% flood regions (swamps and areas affected by floods), and less than 1% mountain vegetation. In recent years, some oil reservoirs have been exploited in Sudan. Our study area is in the oil production fields, located in the south region and lies between latitude 9 to 11 N and longitude 28 to 31 E (see Figure 1). Heglig, Unity, Munga, Bamboo, Diffra, Elnar, Neem, Eltoor, and Toma South are the main oil fields located in Muglad basin, South Kordofan state, and operated by Greater Nile Petroleum Operating Company (GNPOC). It is located within Um Rawaba formations and Muglad basin, which is the largest of the Central African rift basins. The two large hydrocarbon accumulations, the Unity and Heglig fields, have combined recoverable reserves of 250-300 million barrels of oil [39]. The climate of this region is hot with seasonal rainfall, which is affected by the annual shift of the intertropical zone. There are two main seasons in this region. The wet season begins about the end of April and ends about the late of November and the dry season, which lasts for the rest of the year. The highest average temperature is 43 °C in April and the lowest is 33 °C in August. The soil type of the study area ranged from silty sand in the north part to black cotton to the south.

Data and Methods
In this section, we explain the four utilized data sets (i.e., GRACE data, oil wells production record, hydrological models, and Sentinel-1 SAR data) and the implemented methods to study mass change and ground surface subsidence owing to oil and water depletions. According to the time span of the available oil data, which is from January 2003 to September 2012, the GRACE and hydrological data sets were chosen to cover the same period. For the SAR data, we were limited to available Sentinel-1 data, which were freely available, and had a good coverage for our study area. Sentinel-

GRACE Data
The Earth's gravity field and shape change with time owing to various processes. Most of the changes in the Earth are secular and periodic. The largest secular mass changes in the Earth are caused by mantle convection, glacial isostatic adjustment, groundwater level change, as well as plate and intraplate motions [40]. Notable periodic variations are the result of the Earth's rotation (e.g., the pole tide with 14 months' time period), seasonal variations (e.g., water level change) caused by variations in atmospheric and hydrologic conditions, and tidal variations with various periods. All of the above-mentioned phenomena can affect the Earth's gravity field, and it can be measured by space geodesy techniques, for example, GRACE mission. The GRACE data greatly complement some climate change-related aspects such as water resource studies because of their medium-wavelength characteristics. Analyses of the GRACE data allow determination of the temporal changes of the Earth's gravity field to unprecedented accuracy. These data can be used for monitoring various types of mass redistributions in the Earth system, including the ice mass loss, sea level change, groundwater resource variations, and mass variation inside the Earth. Therefore, the GRACE mission is able to provide important information about land hydrology, that is, the sum of land water storage, which consists of surface and groundwater.
The existing global hydrological models suffer from poor quality and lack of in situ data in some areas. One of the goals of the GRACE mission with global coverage was to model the hydrology signals and water mass transport to overcome these shortcomings and improve the existing hydrological models. Different institutes are working on the GRACE data for such purposes, such as the Center for Space Research (CSR) at the University of Texas, Austin; Jet Propulsion Laboratory (JPL) in the USA; and the Deutsches GeoForschungsZentrum, German Research Centre for Geosciences (GFZ) in Potsdam, Germany. The time-varying gravity models are useful for studying mass transport in the Earth. The 30-day series of these time-varying models have been released by the above-mentioned processing centers. Many studies have been conducted to compare the GRACEbased ΔTWS with the hydrological models. Wahr et al. [27] used the output from hydrological, oceanographic, and atmospheric models to estimate the variability in the Earth's gravity field owing to those sources and developed a method for estimating surface mass change from the GRACE solutions. The GRACE can recover changes in the ΔGWS at scales of a few hundred kilometers and larger and at time scales of a few weeks and longer, with accuracies approaching 2 mm in water thickness over land [27]. Yeh et al. [41] studied the TWS change using this mission and concluded that the GRACE-based method of estimating monthly to seasonal GWS changes performs reasonably well at the 200,000 km 2 scale of Illinois. Chen et al. [42] investigated the low degree harmonics of the GRACE solutions and found out that GRACE-based water storage changes are in agreement with estimates from NASA's GLDAS [36] in the Mississippi, Amazon, Ganges, Ob, Zambezi, and Victoria basins.

GRACE Data Analysis and Groundwater Storage Estimation
For this study, the most recent released GRACE monthly solutions (RL06) processed by the CSR (ftp://podaac-ftp.jpl.nasa.gov/allData/grace/L2/CSR/RL06/) were used for quantifying monthly ΔTWS over the period of January 2003 to September 2012, using the time series analysis explained in [43]. Because, using GRACE data, one cannot differentiate between the various surface water storage signals (such as soil moisture, groundwater, and precipitation), the obtained secular rate should be filtered by external hydrological models. The integration of GRACE data with the land surface hydrological models, for example, the GLDAS model, allowed to extract the individual components from the equivalent water storage estimated based on the GRACE solutions and lead to the determination of the ΔGWS.
The GRACE-based ΔTWS was estimated after applying the following steps. Degree one coefficients were added based on the GRACE Technical Note #13b. Degree 2 coefficients were replaced by those introduced in TN-11_C20_SLR_RL06 (GRACE Technical Note TN-11), consistent with the GRACE SDS recommendations. The maximum degree/order of 96/96 was considered, and seasonal signals were not removed. Furthermore, the gravity data collected by the GRACE require smoothing to reduce the effects of errors present in short-wavelength components. Various methods have been proposed to filter the data [33,44]. For example, isotropic Gaussian [27] and non-isotropic [45] filters are the most popular methods. However, none of these methods account for correlated errors in the data.
Owing to the Earth's mass redistribution (e.g., water/oil depletion), a time-dependent change in the gravity field causes changes in the harmonic coefficients, which can be described by the residual spherical harmonic coefficients nm T  that are obtained by subtracting the mean value nm T from the GRACE monthly gravity solutions, where nm T are the numerical values for the potential coefficients. The absolute disturbing potential T and its anomaly ( T  ) can be estimated as follows: where GM is the geocentric gravitational constant (i.e., the product of Newton's gravitational constant and the total mass of the Earth, including the atmosphere), R is the Earth's mean radius, nm Y are the surface spherical harmonics of degree n and order m, and max n is the maximum degree of harmonic expansion. Using linear regression analysis, for the repeated GRACE satellite tracks, secular trend of the disturbing potential can be determined as follows: where nm T   are the secular changes of the residual spherical harmonic coefficients. According to Wahr et al. [27], the GRACE-based secular trend of the ΔTWS can be achieved by where w  is the density of fresh water;  is the Earth's mean density; and n k are load Love numbers of degree n, which can be modeled based on some Earth models [46,47]. For extracting ΔGWS, a hydrological model, for example, GLDAS, should be used to estimate for the ΔSWS (surface water storage anomalies) and ΔSM (soil moisture anomalies). Following the below formula, ΔGWS can be estimated.

GWS TWS SWS SM
where MATLAB scripts were used to analyze GRACE solutions and hydrological models.

Hydrological Data
Land surface model output or hydrological models such as GLDAS, World-Wide Water Resources Assessment (W3RA), and WGHM are very essential tools to compensate for the disturbance factors that contaminate the GRACE-based estimates of ΔGWS (e.g., soil moisture, precipitation). The models provide global information on land surface meteorological and hydrological status (e.g., surface temperature and soil moisture) with temporal resolutions based on integrating satellite and ground-based observations. Different land surface models are currently in use, for example, GLDAS Version-1 [36], specifically Mosaic, NOAH, variable infiltration capacity (VIC), and the common land model (CLM). The GLDAS Version-2 contains four models, NOAH, VIC, CLM, and the catchment land surface model (CLSM). The GLDAS models have different spatial resolutions, for example, 1.0 and 0.25 arc-degrees for NOAH and 1.0 arc-degree for the others. In addition, the layers and depths are varied, for example, NOAH has four layers ranging from 0 to 2.0 meter depths [36]. The GLDAS archived data can be accessed through http://disc.sci.gsfc.nasa.gov/. The WGHM is a 0.5 arc-degree global hydrological model, which was developed at the Centre for Environmental Systems Research of the University of Kassel, Germany, in cooperation with the National Institute of Public Health and the Environment of the Netherlands (RIVM) [37].
In this study, four different data sets of monthly hydrological land surface models, the 1.0 arcdegree of GLDAS data, that is, NOAHv1.0, NOAHv2.1, CLSMv2.0 (Version-1 and -2 respectively), and the 0.5 arc-degree WGHM, were utilized to remove the non-groundwater components from the GRACE-based ΔTWS, over the period of January 2003 to September 2012. Figure 2 shows the secular rate of the sum of ΔSM and ΔSWS obtained from the GLDAS-CLSM model over the study area.

Oil wells Production Data
Sudan started exporting oil in the year 1999, with an official production of about 400,000 barrels per day until mid-2006, and reached about 500,000 barrels per day at the end of 2009 with a proven reserve value of about 5 billion in 2007 [39]. The production and exploration have been limited to the Muglad and Melut basins (see Figure 1a). The Heglig, Unity, Munga, Bamboo, Diffra, Elnar, Neem, Eltoor, and Toma South oil fields are the main ones located in the Muglad Basin, South Kordofan state (see Figure 1b). The data for this study were provided by the Greater Nile Petroleum Operating Company (GNPOC), who operate the oil fields. Oil extraction normally generates large volumes of produced water-approximately 2-10 barrels of water is associated with every barrel of oil produced globally [39]. The oil and water separation process for the nine fields is performed in Heglig Facilities. A daily oil/water production record in m 3 unit, for the nine fields covering the period of Jan 2003 to Sep 2012, was used in this study. Figure 3 shows the monthly total oil and water extractions in the fields. As can be seen in Figure 1b, Heglig has the higher oil production rate, while Toma South and Unity fields are in the second and third places, respectively. Owing to political issues and the separation of South Sudan, oil production dramatically decreased after the year 2012. The estimated associated water to the oil in our fields, for the mentioned period, is about four barrels of water for each one of oil. Thus, as can be seen in Figure 3, almost equal amounts of oil and water are extracted at the beginning of the production, and then gradually, water mass extraction became dominant.

Sentinel-1 Data and Analysis
Synthetic aperture radar (SAR) is an active side looking radar system that has been used in Earth remote sensing and geodesy for more than three decades. Interferometric SAR (InSAR) and differential InSAR (DInSAR) methods use the phase difference between radar images acquired at different times and geometries and allow generation of digital elevation models and measurements of the centimeter-scale earth surface movements [10]. The DInSAR method is mostly used when measuring large deformations are induced, for example, by earthquakes or volcano activities. For small deformation, it has considerable limitations, where the most important one is the atmospheric phase contribution, which can be mitigated by applying the so-called permanent scatterers interferometry (PSI) [38]. PSI exploits multiple SAR images acquired over the same area, and uses appropriate data processing and analysis procedures to separate the displacement phase from the other phase components [48]. The main outcomes of the PSI analysis include the deformation time series and velocity estimated over the analyzed area. The larger the number of available images, the better the quality of the PSI deformation velocity and time-series estimation, for example, at least 15-20 images for C-band sensors [38]. The PSI method is applicable in many areas, such as disaster monitoring and risk assessment, caused as a result of landslides and volcanic eruptions, earthquakes, ice sheet motion and glacier flow, and ground subsidence due to subsurface clay deposits [49] or oil and water depletion [12,50].
In this study, freely available C-band Sentinel-1 A and B data from the European Space Agency (ESA) were used to detect and quantify ground surface deformation related to the groundwater and oil level change for nine oil fields in Muglad Basin, Sudan (Figure 1). For the data processing, SARPROZ software [51] was used, where the single look complex (SLC) images are co-registered to a single master. Then, the SRTM 3-arc second digital elevation model (DEM) and the precise orbits for each image were used to remove the topographic and flat earth components, respectively. Reference points (R) were chosen among the selected persistent scatter candidates (PSCs), which are relatively unaffected by deformation. Permanent scatters (PS) are the strong reflecting objects that are dominant in a cell. Identification of the (PSCs) and PS was performed based on the amplitude stability index and the temporal coherence [52]. Stable pixel selection depended on the use of amplitude dispersion index [53], which can be computed as in Kampes [10]: where A  and A m are the standard deviations and the mean of the amplitude, respectively.
Therefore, a pixel with high reflection has a low A D . The amplitude stability index in the software is defined as . Temporal coherence explains the mean difference between the observed and modeled phase, for each pixel and interferogram [10], and it is a reflection of the standard deviation of the final time series, which is a measure of how well the time series fit its linear regression line.
Three data sets consisting of 30, 30, and 32 single look complex (SLC) images on descending and interferometric wide swath (IW) mode, and single polarization (VV) (see Table 1) were processed and analyzed over the studied oil fields. The images cover both dry and rainy seasons of the periods of Jan 2016-March 2019, Oct 2016-April 2019, and Nov 2015-June 2018, respectively. Figure S1 (see Supplementary Materials) shows the image pair selection for processing of the three data sets. Product type Single look complex (SLC) Polarization VV The first data set, including 30 images, was used to analyze Heglig and Bamboo's oil fields (see Figure 1b), where 0.80 and 0.70 amplitude stability indexes were utilized for the PSCs and PS selection, respectively, and a temporal coherence of 0.60 was used for masking purposes. The second set of data, including 30 images, was used to analyze Diffra (West of the field) and Neem (North of the field) fields (see Figure 1b) where 0.85 and 0.70 amplitude stability indexes and 0.70 temporal coherence were utilized for Diffra's PSCs and PS selection and masking, respectively. Then, 0.75 was chosen as an amplitude stability index and temporal coherence mask for the Neem field. The third set of data, including 32 images, was used to analyze Unity, Munga, Eltoor, Toma South, and Elnar's oil fields (will be referred to as Unity-area), which lie in the South part of the oil fields. The 0.76 and 0.70 amplitude stability indexes and 0.65 temporal coherence were utilized for the PSCs and PS selection and for masking, respectively.

Results and Discussions
In this section, we firstly report our results of GRACE data analysis and the relation between the obtained ΔGWS and the oil production data, and then we continue with PSI results and the estimated land subsidence rate over the selected oil reservoirs owing to oil and water extraction.

GRACE-Detected Mass Changes
We utilized the GRACE monthly solutions processed by the CSR throughout Jan 2003-Sep 2012 to estimate the rate of ΔTWS within the region of interest. Different smoothing isotropic and nonisotropic filters were applied to filter out the effect of the correlated noises present in the monthlynormalized spherical harmonic coefficients. We utilized both DDK filters (non-isotropic) [34,45] and a conventional Gaussian (isotropic filter) spatial smoothing function [35] with 300 and 500 km smoothing radiuses. The estimated changes of ΔTWS over the area are related to the variations in both groundwater/oil storage and other hydrological components, for example, soil moisture. To quantify the ΔGWS variations over the area, the estimated hydrological signal is subtracted from the estimated GRACE-based ΔTWS (see Equation 4). Figures 4a and b show the resulting rate of the ΔTWS change after applying the DDK1 filter to the GRACE solution and the rate of ΔGWS variations based on the GRACE solutions and the GLDAS-CLSM model over the region, respectively. Through the figures, the negative values illustrate mass loss in the study area, while the positive values indicate mass gain.  Figure 4b reveals that the estimated rate of the ΔGWS change over the selected oil fields area is negative (i.e., losing mass). This situation also can be seen in the Northern region, which contains the Nubian Sandstone aquifer, which covers about 2.2 million km 2 of North-East Africa and extended in Sudan, Egypt, Libya, and Chad (see Figure S4 in the Supplementary Materials). CEDARE report [54] showed that there is no significant groundwater recharge for this aquifer and it will be fully depleted in the future, without good planning for extraction. Furthermore, thousands of irrigation pumps located in the East El Owainat, South-West Egypt, which use the nonrenewable Nubian groundwater aquifer, might be the reason for these negative signals, in addition to the artificial Toshka Lake, which was created by the diversion of water from Lake Nasser through a canal into the Sahara Desert. South to Sudan-Egypt border and latitude of 20 degrees, the Nubian aquifer, which extends into Sudan, Chad, and Libya, is in a near-steady state. The dark negative signals in the East (Ethiopia), the South (Uganda), and the South-West, around the Central African Republic and the Democratic Republic of Congo, show significant GWS loss.
One of the goals of this study is to investigate different filters and hydrological models to study mass depletion in the study area in Sudan. As already mentioned, four different hydrological models were used to remove the non-groundwater components from the GRACE-based estimated ΔTWS. Considering the NOAH models, the ΔSM and ΔSWS, which include the surface runoff water, were computed, whereas in the CLSM and WGHM models, the shallow groundwater is also taken into account. According to Vrbka et al. [55], depths to groundwater in Sudan vary from only a few meters in the Nile valley to >250 m in remote and elevated areas. Furthermore, the reported depths for drinking water wells in the Heglig field are around 50 m. Hence, the groundwater in this region may not be counted by the CLSM and WGHM models, as they may normally count for the shallow groundwater regions. The availability of in situ data in the Muglad basin (i.e., oil data) is the reason for performing this study in this area, where Heglig field coordinate (see Figure 1b) was used for the GRACE and hydrological time series data collection (hereafter, it is called Heglig-area). Figure 5 shows the time series based on the four hydrological models and the GRACE monthly solutions using the DDK1 filter in terms of the (ΔSM + ΔSWS) and ΔTWS, respectively, in the Heglig-area. Moreover, the estimated ΔTWS time series based on applying other smoothing and destripping filters to the GRACE solutions are displayed in Figure S2 in the Supplementary Materials. It should be noted that hydrological models differ from each other mainly in terms of assimilation algorithms and input data. For instance, WGHM uses only Global Precipitation Climatology Centre (GPCC) rainfall data as a climate forcing for its simulations. In contrast, GLDAS uses a combined product is derived from terrestrial data and a number of satellite-derived observations [36]. Moreover, to estimate soil moisture variations from the GLDAS, one needs to apply the sum of its four soil moisture layers (0-2 m depth), whereas WGHM represents soil moisture for the complete soil profile in a single layer whose depth refers to the rooting depth of the vegetation cover. Therefore, for instance, for areas classified as bare ground (root depth 0.1 m), the total soil moisture in WGHM could be highly underestimated. Within the GLDAS products, the NOAH model uses the Modified IGBP MODIS 20-category vegetation classification and the soil texture based on the Hybrid STATSGO/FAO datasets. However, the CLSM utilizes the mosaic land cover classification, together with soils, topographic, and other model-specific parameters that were derived in a manner consistent with that of the NASA/GMAO's GEOS-5 climate modeling system [36]. Moreover, the CLSM uses topographically-derived catchment as the land surface element, instead of a grid in traditional land surface models (LSMs). Within the context of comparing NOAHv1.0 with NOAHv2.1, two major issues were found in the GLDASv1.0 forcing fields that have been addressed in GLDASv2.1 product. First, the Air Force Weather Agency's AGRicultural METeorological modeling system (AGRMET) shortwave downward radiation flux was sharp during certain years. Second, there was a dramatic change in precipitation in certain locations starting in 2009. Furthermore, comparisons of GLDASv1.0 radiation and precipitation fields revealed that GLDASv1.0 had high bias relative to the well-validated surface radiation budget (SRB) dataset, and GLDASv1.0 precipitation had low bias relative to the Global Precipitation Climatology Project (GPCP) dataset. Similar biases were observed compared with GLDAS-2.0, whose radiation fields were bias corrected to the SRB dataset (one can find more detailed information in readme document for NASA GLDAS version 2 data products).
The pattern of the four hydrological signals over the region may greatly depend on the climatic status, where they show a bimodal behavior with a rainy season that lasts for about 5-6 months, with a peak in August, and dry season for the rest of the year, with a peak in about February and March. Table 2 provides a comparison between the utilized hydrological models in terms of the correlation coefficients between the estimated (ΔSM + ΔSWS) time series over different periods. This comparison reveals the similarities and differences between the time series estimated based on models, which is highly affected by the lack of and scarce and sparse distribution of in situ hydrological and meteorological measurements in the study region.
A closer inspection of Figure 5 and Table 2 reveals that there is a relatively good agreement  between different time series over 2003-2009, where the NOAHv2.1 model shows a better performance with a correlation of 0.76 with the WGHM model. Then, after 2010, the NOAHv2.1 time series shows a notable deviation from the other models. Moreover, over the whole period 2003-2012, except NOAH 2.1, a better agreement (about 0.52) is shown between CLSM, NOAH1.0, and WGHM models. The peak shown of 2008 may be attributed to the reported high rainfall season that year. In addition, the lower peak pattern in 2009-2010 agrees with the reported low rain season for the same year. Table 2. Correlation coefficients between different hydrological models in terms of the estimated (soil moisture anomalies (ΔSM) + surface water storage anomalies (ΔSWS)) time series over different periods  Table 3 shows the estimated rates of ΔGWS changes using the four hydrological models and different filters. We applied isotropic Gaussian filter, with the half-wavelength of radius of 300 and 500 km, as well as non-isotropic DDK filters to overcome the north-south strips noises. DKK1 is the strongest and DDK8 is the weakest filter. Accordingly, one can see that the estimated uncertainty increases as for applying weaker filters, that is, the weaker the filter, the higher the uncertainty. The results show that using different filters and hydrological models will result in different patterns of rates of ΔGWS changes. For instance, using different filters, the NOAHv2.1 hydrological model is the only one among the others that results in a positive rate of ΔGWS change. Moreover, the CLSM model provides a higher negative rate of ΔGWS change when compared with others. Table 3. The rate and their uncertainties of groundwater storage anomaly (ΔGWS) change based on the Gravity Recovery and Climate Experiment (GRACE) products and hydrological models for different filters. Unit: cm/yr (Uncertainties are computed based on the 95% confidence level). To evaluate different filters and hydrological models, one can rely upon the available in situ data (i.e., oil production), as the total mass extraction should agree with the estimated ΔGWS. Table 4 summarizes the correlation coefficients between the monthly total oil extraction and the estimated ΔGWS change using different smoothing filters and hydrological models in the Heglig-area oil fields. We used the well-known Pearson's equation [56] to quantify the similarity between in situ-achieved oil and water extraction and the GRACE-based ΔGWS estimates. The estimated ΔGWS based on the DDK1 and CLSM model shows a good agreement with the local mass extraction (see Figure 3), that is, oil and water depletion (with a correlation coefficient of about 0.74).  Figure 6 illustrates the estimated ΔGWS based on the different hydrological models and applying the DDK1 filter to the GRACE monthly solutions in the Heglig-area. Similar to Figure 5, the obtained ΔGWS time series from different hydrological models are generally agreed between 2003 and 2009, however, after 2009, the estimated ΔGWS using the NOAHv2.1 model shows a significant deviation from the other models.

PSI-Detected Land Subsidence
The availability of the proper SAR satellite data over the given region and revisiting time affects the capability to estimate deformation phenomena over time [57]. Despite the suitability of L-band images for vegetated areas in terms of penetration depths, there were not enough free L-band images to cover the study region for such a PSI analysis. As such, only Sentinel-1 data were used for this study. As can be seen in the graph ( Figure S1), there are a lack of data during periods of 23 Jan 2016 and 21 July 2016 (six months) for the first data set, in addition to some sparse periods (2-3 months) between the acquisitions of the three data sets despite the short revisit time of Sentinel-1 sensors, which is 6-12 days.
The three SAR data sets were analyzed over the nine oil fields, where the PSI technique was applied to detect relative ground deformation due to oil and water depletion, using the processing criteria in Section 2.4. The line-of-sight (LOS) displacement rates at the selected PS points in Heglig and Bamboo (Figure 7), Neem ( Figure S6), Diffra (Figure 8), and the Unity-area fields ( Figure S7b) were computed relative to their reference PS point (cf. Figure 7).   Table 5 summarizes the cumulative displacement LOS, temporal coherence, velocity, and velocity standard deviation of the selected PS points (i.e., those have maximum velocity), including their references for Heglig and Bamboo, Neem, Diffra, and Unity-area oil fields, respectively. Furthermore, Figures 9 and Figure S7 illustrate the cumulative LOS displacement in millimetre for the selected PS points relative to their corresponding reference points. Land covers of Neem oil field, which is located in the northern part of the study area, are mostly silty sand and less vegetated with limited surface structure (i.e., most of the PS points are laid on the ground) (see Figure S6). Reference point, R (N-568), at the top of the figure, can be seen among the very stable PS points (i.e., green color), and shows a stable and consistent time series ( Figure S6). The time series of PS points N-586, N-653, N-669, N-741, and N-761 shows a continuous subsiding rate with notable deviations in May, June, and July that can be attributed to the rainy season.
To the south of Neem, where Heglig, Bamboo, and Differa fields are located, the land is more vegetated and soil types are more clay. Heglig is the head quarter of the nine fields and it contains most of the oil infrastructure (oil central processing facilities (CPFs), airport, and the main camps). This is why the density of the PS points is larger compared with that of other fields (see Figure 7). The reference point R(H-124) in Figure 7 can be seen among the very stable PS points (i.e., green color), and also shows stable time series, as can be seen in Figure 9a. In this figure, the time series of point (H-178) illustrates a consistent subsidence rate of 10.53 mm/yr. The time series of the PS point (Heg-166) with an estimated 24.47 mm/yr subsidence rate shows notable deviations during months of April and May, which can be attributed to the rainy season. The three PS points (Bam-6, Bam-8, and Bam-9) are in Bamboo field, and show consistent time series generally with small changes in the rainy seasons. Therefore, the reasons for such deformation in this area can be attributed mainly to the oil extraction and the soil type. The use of pile foundation type in building constructions, which is dominant in Heglig field, can greatly reduce the soil type effect. Thus, the main reason for such displacement in Heglig, which has a higher production rate, probably could be assigned to oil and water depletion. The reference PS point in Diffra field (Dif-81) shows more temporal varieties in May, June, and October (i.e., the rainy season) (Figure 9b).
To the south, in the Unity-area fields, the most vegetated area and large swamps limited our PSI analysis. Figure S7b shows the reference PS point, R (Uni-40), which is less stable compared with other references. For these fields, the effect of the soil type and swamps, and the lack of stable reflectors (e.g. buildings) on the PS time series analysis are visible.
The comparison of the subsidence rates for different oil fields (Table 5) clearly shows clearly the Heglig oil field is the most deformed one, with an accumulative displacement of about 77.6 mm during a three-year period. Its higher production rate may be the reason for such large subsidence, in addition to the dense PS points, which allow for closer inspection. Moreover, and according to Castelazzi [58], InSAR applications to groundwater depletion are limited to reservoir systems sensitive to measurable deformation (clays and silts), which is the case in the study area.

Potential and Limitations for Linking GRACE and PSI Results
Despite their different resolution and application, GRACE and SAR can provide complementary data for monitoring of water and oil reservoirs [58]. The water and oil depletion causes changes in mass distribution, and the consequent reservoir compaction causes ground surface subsidence, which can be detected by InSAR (e.g., PSI analysis in our case). According to [58,59], the generated PSI-based groundwater depletion map can be used as downscaling data to partially mitigate GRACE resolution limitation. For such a process, additional in situ data are required (e.g., the vertical and horizontal variability of sediment compressibility) to achieve such an inversion from PSI-based displacement map to volume of depleted GWS. However, this analysis was not included in this study owing to a lack of such data.
In the context of this study, and owing to the lack of SAR data for our study area before 2015, we could not have the temporal overlap with other data sets, that is, GRACE and oil wells data.
However, owing to the continuous oil production between 2002 and 2019 ( Figure S5), the analyzed trend of subsidence over oil fields resulting from our PSI analysis can be projected approximately for different periods in which GRACE and oil wells data sets are available. However, with the accumulation of Sentinel-1 or other SAR data together with the new data collected by, for example, GRACE-FO, a more detailed analysis and link between results will be possible.

Conclusions
In this paper, we carried out a comprehensive study about the relation of oil production, GRACE-based ΔGWS estimation, and monitoring ground subsidence in response to extraction cycles in nine oil fields in Sudan. We used GRACE monthly solutions (RL06) processed by the CSR processing centre to assess the oil and water depletion in the selected oil fields. In this study, we employed both non-isotropic (DDK filters) and isotropic (Gaussian) filters to reduce the effect of the correlated noises in GRACE data. We computed the rates of ΔGWS changes (due to the oil production) by utilizing different hydrological models, namely, GLDAS-NOAHv1.0, GLDAS-NOAHv2.1, GLDAS-CLSMv2.0, and WGHM models.
Our results show that the achieved rates of ΔGWS changes are highly dependent on the choice of decorrelation filter and the hydrological models. Our analysis based on applying DDK1 decorrelation filter and utilizing the GLDAS-CLSMv2.0 model results in a high correlation (0.74 correlation coefficients) between the estimated ΔGWS and oil and water extraction in the Heglig-area oil fields. Except for the resolution of GRACE data, the main source of uncertainty in the estimated ΔGWS is related to the hydrological models, which reveal a strong relation with the rainfall trend in the region. According to the achieved results, GRACE data can help to monitor the mass changes owing to heavy oil productions, if we could solve remarkable uncertainty in hydrological models in the future.
We also computed the rate of land subsidence owing to the oil extraction in different oil fields using SAR data and the persistent scatterer interferometry (PSI) technique. For instance, we observed 77.6, 60.1, 35.8, and 15.4 mm subsidence in Heglig, Neem, Diffra, and Unity-area oil fields during a three-year period (i.e., Nov 2015-Apr 2019), which correspond to 24.47 ± 0.85, 23.8 ± 0.96, 14.2 ± 0.85, and 6 ± 0.88 mm/yr, respectively. This surface subsidence can be assigned to the pore pressure drop owing to the vertical compaction of the reservoir.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/12/11/1792/s1, Figure S1: Image graph for the 30, 32, and 30 descending sets of images, (a) for Diffra and Neem's sites; (b) for Unity, Munga, Eltoor, Toma South, and Elnar's sites; and (c) for Heglig and Bamboo's sites, respectively. The figure shows acquisition date (horizontal axis), chosen master, and the perpendicular baselines in meter unit; Figure S2: Computed TWS based on the different hydrological models and GRACE using different fillters in Heglig area; Figure S3: GWS obtained from subtraction of GRACE (different filters) and different hydrological data; Figure S4: Groundwater resources in Sudan and South Sudan; Figure S5: Oil production in Sudan and South Sudan; Figure S6: Displacement trend of the descending PS points at Neem oil field relative to reference point R (N-568). Unit: mm/yr. Period: Oct 2016-Apr 2019; Figure S7: Cumulative displacement in (mm) of the selected PSI points relative to reference point R in (a) Neem and (b) Unity area oil fields.
Author Contributions: N.A.A.G. wrote the manuscript and processed and analyzed the data. H.A. carried out data processing and helped to draft the manuscript. M.B. designed the study and helped to write and prepare the manuscript. He also contributed to data processing. F.N. contributed to data processing and edited the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.