Using InSAR Time Series to Monitor Surface Fractures and Fissures in the Al-Yutamah Valley, Western Arabia

Abstract

Western Arabia routinely experiences geophysical phenomena that deform the surface of the earth in a variety of ways. These phenomena include earthquakes, volcanic eruptions, sinkholes, and earth fissuring and fracturing. We perform a time-series analysis of interferometric synthetic aperture radar (InSAR) observations derived from the ESA Sentinel-1 radar satellite constellation to map regional surface displacements in western Arabia as a function of time. We rely on InSAR products generated by the JPL-Caltech ARIA project to detect regions with short wavelength anomalies, and then manually reprocess InSAR products at a higher resolution for these regions to maximize spatial and temporal coverage. We post-process InSAR products using MintPy workflows to develop the InSAR time series. We report short wavelength anomalies localized within alluvial valleys across western Arabia and find a 5 cm/year line-of-sight surface displacement within the Al-Yutamah Valley. Part of the observed subsidence is correlated with surface fractures that developed in conjunction with severe rainfall events in regions characterized mainly by alluvial sediments at the surface. Regions of observed subsidence that are not associated with any surface fractures or fissures are correlated with the presence of basalt layers at the surface. Both regions are subject to groundwater exploitation. The observed subsidence is inferred to be driven by groundwater withdrawal perhaps modulated by the presence of a preexisting depositional environment (e.g., paleo-lake deposits) that promotes unconsolidated soil compaction.

1. Introduction

The western Arabian margin (Figure 1A) has experienced a rich history of volcanism during the last 20 Ma [1]. However, these volcanic events are not the only potential geologic hazard that affects this region. Surface fractures and sinkholes within the Arabian margin are also rapidly becoming a challenging hazard. These hazards are frequently linked to surface subsidence. Thus, identifying such regions of subsidence can provide an indication of potential future hazards [2,3,4,5,6]. Recent studies of surface deformation within Arabia investigate surface fissures and fractures, explore their underlying mechanism, and address mitigation of their potential risks [7,8,9,10,11,12,13,14,15,16,17,18,19,20]. These studies used multiple methods, including structural geology, soil analysis, and interferometric synthetic aperture radar (InSAR). Among these approaches, InSAR is becoming one of the most powerful tools that can be used to monitor ground movement over time. This approach relies on determining phase differences in repeat synthetic aperture radar (SAR) images to quantify spatial variations in surface displacement (e.g., [21]). The InSAR phase is a superposition of various signals, including all geophysical contributions to ground displacements and propagation delays accrued in the atmosphere. Thus, we must pay particular attention to separating the various contributions to the measured phase change in order to isolate the geophysical signal of interest.

In this work, we apply an InSAR time series analysis of Sentinel-1 data to cover wide areas (250 × 1000 km) of western Arabia at a 90 m spatial resolution and with short revisit times (6–12 days on average). This regional survey allows us to detect the short wavelength anomalies that correspond to local ongoing deformation within the Arabian margin. Then, we reprocess the Sentinel-1 data at a much higher (15 m) spatial resolution for a limited region where we have identified short wavelength deformation anomalies. We use this local analysis to optimize spatial and temporal resolution for the Al-Yutamah Valley in comparison with our regional survey, to investigate domains of earth fissuring and fracturing, and to monitor their activities over time. We develop time series, velocity maps, and fracture event spanning maps. Here, we present both the regional and local analyses, correlate the observed deformation fields with surface fractures and fissures, discuss the possible underlying mechanisms, and forecast future hazards.

Study Area

Within western Arabia lies an exposed Neoproterozoic block (Figure 1A) known as the Arabian Shield, with a crystalline Precambrian basement, Cenozoic alkalic volcanic fields (Harrat), and small sedimentary basins and grabens (e.g., [22]). The Al-Yutamah Valley (Figure 1B) is located in western Arabia (23.8852 N, 39.7524 E) within the Arabian shield, in the vicinity of Harrat Rashid [23] or Al-Madinah [24]; which is the northern part of the Harrat Rahat volcanic field, about 100 km from the Red Sea coast. The study area (Figure 1B) is dominated by basaltic lava flows of the middle–Pleistocene age that bound the area from the eastern and the northern sides. Precambrian basement is exposed in the western and southern regions, while surficial deposits (Alluvium—Quaternary) present as mud/loam flat around the periphery of volcanic rocks [25]. The Al-Yutamah Valley was formed by E–W Cenozoic dextral strike slip faulting (Figure 1B) [26], which has been suggested to be associated with opening of the Red Sea [27]. The other major structural features of the area generally trend in NW–SE direction, consistent with the regional structural direction and believed to be part of the Najd strike-slip fault system (NFS) [26]. However, the pattern of small-scale faulting changes strike in the Al-Yutamah basin to a N–S direction, as seen in the field within the Quaternary deposits [8].

2. Methodology

2.1. Geologic Map

We have considered various geological maps for the study area from the Saudi Geological Survey [30,31] and the United State Geological Survey [25,32]. While the most recent detailed map published by USGS [25] has the highest resolution boundary between the geological units, our area of interest has not been completely mapped (about 30% of the area has been covered). Thus, these resources were supplemented with our field observations and high-resolution satellite images (Figure 1C) to create a complete map for the study area (Figure 1B). The geological map of the Al-Yutamah contains three main units (Precambrian shield, Pleistocene basalts, and Quaternary deposits) with the basalts unit further subdivided into three units based on the recent observation and classification of [25].

2.2. Weather Model and Ancillary Weather Data

To correct the InSAR time series for tropospheric delay effects, we rely on the weather models provided by Copernicus Climate Change Service (ERA5) [33]. The ERA5 is the fifth generation of the European Centre for Medium Range Weather Forecasting (ECMWF) [34], an hourly analysis for the global climate and weather. This model has been regridded on a regular lat–lon grid of 0.25 degrees. We also compare our InSAR time series and rainfall in the Al-Yutamah Valley (See Section 3.2.1). The rainfall data are taken from World Weather Online (https://www.worldweatheronline.com/aboutus.aspx; accessed on 20 February 2022, last accessed on 31 March 2022).

2.3. InSAR Data and Method

We use radar imagery acquired by the Copernicus Program Sentinel-1A and 1B (track 14) synthetic aperture radar (SAR) satellites developed and operated by the European Space Agency. We consider InSAR images at both regional and local scales. We start with a regional scale (100 s of km) to detect short wavelength anomalies that correspond to local ongoing deformation and consider local scales (10 s of km) to maximize spatial and temporal resolution for local regions of interest. All our InSAR processing relies on the open-source InSAR Scientific Computing Environment (ISCE; [35]) software.

At the regional scale, we use the so-called “standard product” produced by the Advanced Rapid Imaging and Analysis (ARIA) system [36]. The ARIA system uses algorithms from ISCE to provide an operational archive of standard Sentinel-1 Geocoded UNWrapped interferogram (GUNW) products, enabling rapid ingestion into post-processing algorithms. These InSAR-derived surface displacement products are reduced from the original resolution down to a 90 m spatial resolution to increase the signal-to-noise ratio. Interferometric pairs are generated between the nearest two acquisitions, leading to temporal baselines of 12 and 24 days. We also produce additional annual interferogram pairs spanning October and May of each year. We use the ARIA-tools open-source package [37] to download and stitch 5102 GUNW products into 252 continuous interferograms from 133 acquisitions spanning the nearly 5-and-a-half-year period between February 2016 and July 2021. The region covered extends from 20N to 28N and 37E to 41E. All our time series are referenced to the date 5 February 2018. We consider the root-mean-square (RMS) residual phase as a measure for our data quality (Figure S1). The original number of acquisitions is 183, where we removed 50 acquisitions using a 24 mm residual threshold. We evaluate the quality of our data by calculating the history of interferometric coherence (Figure S2, left panel). The average coherence is high, with values generally above 0.85. We developed a redundant network for western Arabia (Figure S2, right panel) with each acquisition connected to an average of three other acquisitions.

We applied four main corrections to the InSAR data:

• Model-based removal of the deformation signal caused by the solid earth tides (SET) [38,39] (Figure S3). There are both significant spatial and temporal variations in the impacts of SET.
• Model-based removal of radar propagation effects associated with temporal variability in the quasi-stratified component of the troposphere [40].
• Estimation and removal of error associated with the inaccurate digital elevation models (DEMs), as the sensitivity to this error is proportional to the perpendicular component of the interferometric spatial baseline (so-called B_perp) and is easily estimated in the time series analysis [41].
• Empirical removal of residual long wavelength (100-km-scale) apparent deformation signals [42]. This removal is recommended for short spatial wavelength deformation signals, which estimate and then remove linear ramps from the displacement time-series at each acquisition on the reliable pixels. In contrast, for regional tectonic deformation signals, this approach is not recommended.

We use the ARIA tools package suite to prepare interferograms for time series analysis in the Miami INsar Time series software in Python language (MintPy) [42]. This workflow has been successfully demonstrated for mapping ground deformation in subsidence and tectonic active regions along the Eastern US [37], US west coast [43], and the Afar region of Ethiopia [44].

We also generate a higher-resolution InSAR dataset, with ~16 m spatial resolution, to further investigate the local deformation around the Al-Yutamah Valley identified from the ARIA products. We use Sentinel-1 ascending track 14 data acquired in Interferometric Wide (IW) mode from 6 October 2014 to 19 July 2021 (202 acquisitions in total). We use the stack processor [45] within the ISCE software [35] for interferogram processing. The stack of SAR images is coregistered using pure geometry with precise orbits and a digital elevation model (DEM). We remove the topographic phase using the SRTM DEM (SRTMGL1, 1 arc second with voids-filled, ~30 m) [46]. We pair each SAR image with its nearest 3 neighboring images in time domain to form 596 small temporal baseline interferograms, multi-look each interferogram by 3 in range, and filter the phase using the Goldstein filter with a strength of 0.2. We excluded 16 interferograms due to low coherence. The interferograms are phase unwrapped using the minimum cost flow method [47]. Increasing the spatial resolution comes at a cost of lowering coherence. We fit a simple temporal model that includes linear and periodic (annual and semi-annual) terms, and consider RMS residual phase in mm as a measure for our data quality (Figure S10). The original number of acquisitions is 264. We identified and removed 62 particularly noisy acquisitions using a 5 mm residual threshold. The coherence history (Figure S11, left panel) of our region shows a lower coherence in comparison with our regional scale results. This decrease is caused by increasing spatial resolution. Nevertheless, the temporal coherence of the Al-Yutumah valley is high (S12), indicating minimum unwrapping error. Moreover, as phase unwrapping error can be detected using the number of interferogram triplets with non-zero integer ambiguity of the closure phase [42], our results (S13) show two triplet integer peaks suggesting 2 out of 580 interferograms have unwrapping errors in the coherent region, which is a minimal amount, and should has negligible effect on our time series estimation. We use the MintPy software [42] to form the InSAR time series from the stack of interferograms using the small baseline approach [48]. We exclude low-coherence interferograms using the coherence-based network modification with a threshold of 0.7. We remove the stratified tropospheric delay using the ERA5 weather reanalysis dataset [33] with the PyAPS software [40], the topographic residual caused by the DEM error [41], and a linear phase ramp for each acquisition to remove the long spatial wavelength component. We use the temporal coherence [49] to eliminate the unreliable pixels with a threshold of 0.7. The noisy acquisitions are identified based on the residual phase root mean squares with a predefined cutoff value (3 times the median absolute deviations) and excluded during the estimation of topographic residual and average velocity [42].

2.4. Ground Truth Data

We conducted multiple field trips in 2021 and 2022 to the Al-Yutamah Valley to map the developed fractures and correlate them with observed InSAR subsidence. We surveyed the fractures that developed in 2018 based on local witnesses and satellite images that were taken by the Worldview 2.0 optical satellite five months after the fracture event. We surveyed the important infrastructure that intersects with our observed InSAR subsidence and visited a couple of sites damaged by the developing fractures.

3. Results

3.1. Regional Analysis of Western Arabia

We discuss here the impact of each of the corrections we considered for the InSAR time series. First, we corrected our phase time series for solid earth tides by removing the calculated solid earth tide displacement in the LOS from our observation. The resultant map (Figure S4) shows regional long wavelength anomalies as the displacement increases gradually westward. However, short wavelength anomalies appear within western Arabia, distributed away from the Red Sea coast, within the plate interior. Second, we remove the tropospheric phase residual, and the resultant displacement map of this correction shows the observed displacement field corrected for the tropospheric component (Figure S5). Comparing this displacement field with the displacement field corrected for solid earth tides (Figure S4) suggests minor changes, as the long wavelength anomalies still sweep across western Arabia. Third, removing the delay corresponding to DEM error resulted in a displacement field (Figure S6) that shows very minor effects on the observed regional long wavelength signal.

A long wavelength phase residual still exists, which may be caused by ionosphere effects, or inaccuracies in the assumed orbits. We remove this long wavelength signal empirically by fitting and removing a bilinear ramp from the phase [42]. We also estimate annual (Figure S7) and semi-annual (Figure S8) components from our displacement time series. Our annual and semi-annual magnitude maps show potential long wavelength signals between latitudes 22N and 26N. These signals are spatially correlated with the observed phase residual magnitude transition from east to west. These observed signals are not validated nor well understood and are thus simply removed from our velocity since we are primarily concerned with short length scale displacements.

The resultant velocity map (Figure 2A), which excludes the annual and semi-annual signals, shows long wavelength velocity variations with small magnitudes. On the other hand, this velocity map also shows many prominent short wavelength velocity variations with larger magnitudes. The uncertainty of our results increases proportionally with distance away from the chosen spatial reference point as shown in the standard deviation of our velocity map (Figure S9). We show three examples of observed subsidence within western Arabia, the Alyutamah, Dhaa, and Shajwah Valleys (Figure 2B–D respectively). The observed subsidence within all the reported examples exceeds 2 cm/y, which is much larger than the noise threshold (Figure S9). Our focus in this paper is on the Al-Yutamah Valley (Figure 2B). For this region, we manually reprocess the InSAR data to increase the resolution and enhance the spatial and temporal coverage.

3.2. Al-Yutamah Valley

3.2.1. Velocity Map

We describe sustained subsidence in the Al-Yutamah Valley (Figure 2B) that has been observed using the regional, low resolution InSAR time series. We use 580 interferograms to develop the displacement time series for the Al-Yutamah Valley. Since we are working at a local scale, we have corrected our phase time series for troposphere and DEM error only. The resultant displacement time series has been used to estimate the annual and semi-annual signals, which are in our case lower than 0.4 cm magnitude. We use a linear temporal model to produce the velocity map of our area. The resultant velocity map (Figure 3, left panel) shows similar subsidence to what we already observed with our previous analysis (Figure 2B), but with much higher resolution and larger magnitudes. The left panel of Figure 3 shows the subsidence in our region, where we mask the lower temporal coherence domains as white. We present the displacement time series within the middle of the subsidence (Figure 3, right panel), which shows significant accumulated linear subsidence with 40 cm over the last seven years. A few minor offsets appear in the displacement time series, they can be addressed by removing the linear trend (analysis to follow).

3.2.2. Fractures

We observed two dominant temporal behaviors in our displacement time series, the nearly constant accumulation of subsidence and discrete offsets that occurred on November 2018. We describe these two components in more detail below. We visualize the discrete offset with our high-resolution dataset by fitting a step function at the beginning of November 2018. The resultant event spanning map (Figure 4, left panel) has been spatially smoothed using a Gaussian filter (half-width 3.0 pixels). We find that both uplift and subsidence occurred within a week, with a peak amplitude of 3 cm, which is much larger than the calculated error (Figure 4, right panel).

We describe four local regions (Figure 5A–D) and correlate them with rainfall accumulations. The time series presented in Figure 5A–D exhibits strong correlation with the rainfall events which occurred during November 2018. Beside each time series, we present the corresponding image that taken in April 2019 by the Worldview 2.0 optical satellite, which is four months after the fracture development/reactivation event. We removed the linear component (detrended) out of our time series in Figure 5A,B,D, while C has no prior linear subsidence. We validated our InSAR time series with field observations.

-

(Figure 5A) Subsidence within the Al-Yutamah Valley is associated with distinct surface fractures. In the right panel, we show the displacement time series for the location indicated by the black triangle in the upper left map to better visualize these offsets. The time series shows accumulated subsidence over the last seven years, with multiple offsets, one of which occurred during November 2018 associated with the fracture that can be seen in the optical image (Figure 5A).

-

(Figure 5B) The associated displacement time series (at the location indicated by the black triangle) has multiple temporal offsets. In the left panel is an image which shows the fracture in green, as it accommodated large volumes of rain and became a suitable environment conducive to plant growth. In the right side of the left panel is a picture taken by one local resident of the fracture zone upon its initiation. The picture exemplifies the size of the fracture and the potential hazard posed by such features. This region subjected to subsidence, albeit at a lower rate than the region presented in Figure 5A, also developed surface fractures in response to the severe rain season in November 2018.

-

(Figure 5C) A region outside of the original subsidence area showing little or no subsidence activities prior to the fracture development that occurred in November 2018. In the left image is the region of subsidence, which revealed a potential pattern which appeared in the form of fractures that connected in circles, which may lead into developing future sinkholes.

-

(Figure 5D) A region experiencing a long period displacement variation with a unique pattern of surface displacements that may reflect subsurface processes. The observed pattern seems to reflect a long-term subsurface process within the valley.

3.2.3. Spatial Gradient of the LOS Displacement Field

We calculate the first spatial derivative of the observed InSAR LOS velocity field (Figure 6) to quantify the distribution of the local deformation pattern due to the accumulative subsidence. We produce this map by computing the gradient amplitude as independent from the LOS direction. The resulting map shows high gradients along the edges of the subsidence area. We present three cross-sections, the first (Figure 6A-A’) shows the gradient amplitude of the LOS velocity field along the Haramain High Speed Railway that connects Jeddah to Madinah. The second transect (Figure 6B-B’) shows the gradient amplitude of the LOS velocity field along the national east–west pipeline that connects Beqaiq to Yanbu. The third transect (Figure 6C-C’) shows apparent gradient amplitude of the LOS velocity field within the valley, where the cumulative apparent gradient amplitude maximum value reaches 15 cm for the period between 2014 and 2021. The zones of high gradient are highly correlated with reported surface fractures (Figure 5). In contrast, although the apparent gradient amplitude of the LOS velocity field exhibits similar magnitude in the adjacent region (along the pipeline (Figure 6B-B’ transect)), there is currently no evidence for surface fractures.

4. Discussion

Using InSAR time series and ground truth observations, we investigated the occurrence and recent development of local fractures within the Al-Yutamah Valley, western Arabia. These hazards are associated with land deformation/subsidence as mapped using InSAR and identify potential areas at risk. We report three main observations: continuous subsidence since 2014 (Figure 3 right panel), development of large surface fracture events (Figure 5A–D), and long period oscillations in surface displacements (Figure 5D).

The long-term linear subsidence pattern (Figure 3, right panel) within the Al-Yutamah Valley suggests a long-lived subsurface process. The observed subsidence distribution prompts many questions about the nature of its causes. The most obvious of which is the role of groundwater exploitation that may also drive compaction processes within a subsurface aquifer [8] (Figure 3, right panel). Multiple studies of different examples around the world [16,50,51,52] show a strong relationship between groundwater exploitation and surface subsidence. Although the subsidence area appeared in the alluvial sediments and basalt flow, the fractures are observed only in the alluvial sediments (Figure 7). The observed pattern (Figure 5D) displays a sinkhole-like shape within the basalt area linked with a long meandering region of subsidence. Prior studies have attributed this subsidence to withdrawal of groundwater [7,8] based on an analysis collected of soil samples from the Al-Yutamah Valley. In the Najran Valley [12] geophysical surveys were conducted across observed surface fissures and they were attributed to water pumping and consequent continuous decline in groundwater, which led to the development of subsurface topography characterized with steep slopes (buried cliff) that are correlated with the observed surface fissures.

Reference [8] suggested another possible process based on their soil analysis for surface fractures in the Al-Yutamah Valley as developed due to the hydrocompaction of soil after severe rainy seasons. Though the local accuracy of rainfall accumulations could be missing due to the lack of weather stations in our region as the rainfall data are extracted from gridded dataset. We find that severe rainfall events are correlated with the observed offsets in our displacement time series. One of these events occurred in November of 2018 (Figure 4, left panel). During that event, the displacement field shows uplift along the regions that results in long term subsidence and subsidence along its boundaries. Furthermore, all the developed fractures are located along the boundaries of the original subsidence (Figure 7). Despite the potential effect of rainfall, part of the subsidence (Figure 5C) does not show any correlation with the repeated periods of severe rainfall, suggesting the possibility subsurface cavities may develop because of excessive withdrawal of groundwater from a subsurface aquifer. At the same time, we observed another pattern downstream at the western end of the subsidence region (Figure 5D) that reveals a long-term uplift/subsidence cycle, which may suggest aquifer related processes—but at this time, we cannot explain the actual process.

An alternative hypothesis relates the subsidence to subsurface rheological changes that promote compaction due to groundwater exploitation, leading to the development of fracture zones and sinkholes [4,6]. This hypothesis is motivated by the fact that the subsidence region has sharp boundaries that suggest sharp contrasts between subsurface regions with different rheological behavior, perhaps a consequence of a paleo-lake/paleo-river depositional environment; this implies that the subsidence region is characterized by mud deposits, which is different from the surrounding sandy deposits. Such a process has been reported elsewhere [53,54,55], particularly within the regions that are subjected to groundwater exploitation.

Regardless of the cause of Al-Yutamah subsidence, the area (Figure 7) is clearly subject to natural hazards that can cause severe damage to the residents and/or the major infrastructures such as the national east–west pipeline and the Haramain High Speed Railway. The cumulative gradient amplitude of the LOS velocity field that is observed in the Valley, where most of the fractures were reported, is similar to what is observed (Figure 7) in the vicinity of the national pipeline. That fractures have not yet manifested themselves at the surface does not mean that we should not expect them in the future if the subsidence continues. In fact, given the accumulation of strain, we may expect larger surface fractures or fissures than what we have already observed in the adjacent valley.

5. Conclusions

We mapped the regional and local velocity fields within western Arabia using a time-series analysis of InSAR observations derived from the ESA Sentinel-1 radar satellite constellation. We used InSAR products generated by the JPL-Caltech ARIA project for regional analysis, and we manually reprocessed InSAR products with 580 interferograms at a higher resolution for local analysis to maximize spatial and temporal coverage. We reported three regions where we observed continued subsidence within western Arabia (Al-Yutumah, Shajwah, and Dhaa valleys) (Figure 2). We focused on the Al-Yutumah valley where multiple fissures and fractures were reported. We applied multiple corrections on the regional and local scales. For the regional, we removed the delay within our time series due solid earth tides, troposphere, and DEM error. Additional long wavelength removal applied to highlight the short wavelength signals in our regional velocity map. For the local scale, we removed the delay corresponding to the troposphere and DEM error only to preserve any potential signal.

We observed continued subsidence within the Al-Yutumah valley with more than 40 cm accumulated displacement (Figure 3). We mapped the spanning event in November 2018 to show the spatial distribution of observed displacement in consequence to rainfall event (Figure 4). We presented the time series for the regions where fissures and fractures were observed. We manifested the strong correlation between the main rainfall event in November 2018 and the observed offset in our time series (Figure 5). We calculated the spatial gradient amplitude of the observed LOS InSAR velocity field to highlight local deformation patterns (Figure 6). The spatial gradient amplitude reaches its maximum value along the boundaries of the subsidence, which shows strong correlation with the distribution of fractures and fissures in the Al-Yutumah valley. We plotted the spatial gradient along three main cross-sections, first along the railway, second along the national pipeline, and third across the valley where we observed surface fissures and fractures (Figure 6). We manifested that the magnitude of spatial gradient within the valley is close to the magnitude along the pipeline, suggesting potential fractures and risk.