Use of Geophysical and Radar Interferometric Techniques to Monitor Land Deformation Associated with the Jazan Salt Diapir, Jazan city, Saudi Arabia

Using integrated Interferometric Synthetic Aperture Radar (InSAR) datasets (Envisat: 2003–2009; Sentinel-1: 2014–2018), local gravity surveys, and passive seismic data, we investigated the environmental hazards associated with the rise of the Miocene Jazan salt diapir (JZD; ~ 2 km2) within Jazan city, Saudi Arabia, and identified areas at risk in its immediate surroundings. Our findings include (1) the JZD outcrop and its northern, southern and western bordering areas have been undergoing substantial uplift (up to 4.7 mm/yr), whereas the sabkhas to the east are witnessing subsidence (up to − 7.5 mm/yr); (2) a low Bouguer anomaly (7.5 mGal) was observed over the JZD relative to its surroundings (8.5–12 mGal) with the steepest gradient along its eastern side; (3) strong and clear horizontal/vertical (H/V) spectral ratio peak and high frequency (5–10 Hz) over the JZD outcrop and areas proximal to its western margin, but areas to the east have a weak H/V peak and low frequency (1.5-3 Hz); (4) drilling confirmed presence of a shallow (4 m) salt bedrock layer west of the JZD and the absence of this layer to its east (up to depths of 60 m); (5) uplift patterns along the diapir margins are indicative of near-vertical contact along the JZD eastern margin and less steep contacts along the remaining margins; and (6) additional near-surface diapirs could potentially be identified in the vicinity of the JZD using our integrated approach.


Introduction
Salt diapirs have been reported in more than 35 basins worldwide Jackson 2006, 2007). The salt diapirs often contain other interbedded evaporitic (e.g., anhydrite or gypsum) and non-evaporitic units. Large salt bodies deposit in four main settings: cratonic basins, synrift basins, postrift passive margins, and continental collision zones (Hudec and Jackson 2007). Salt diapirs develop when the less dense and incompressible salt layers are buried under the more compact and dense overburden sediments. This process creates a density inversion, and over time the less dense and buoyant salt layers may pierce through the overlying rocks. The ascent of the diapir is also driven by differential loading of the overburden and regional tectonics (Chapman 1983;Harding and Huuse 2015;Jackson and Hudec 2017;El Rabia et al. 2018). For buoyancy alone to cause the upward migration of the diapir, the clastic overburden has to be buried to at least 1600 m, but more commonly to 3000 m (Jackson and Hudec 2017). The diapirs are usually, but not always, associated with a caprock (Goldman 1933). Caprock forms during the ascent of the diapir, when some of the halite dissolves and gets chemically altered leaving behind a hard competent layer largely made up of anhydrite or gypsum that protects the underlying weaker salt layer (Goldman 1933;Jackson and Hudec 2017). The caprock is often unevenly distributed, especially once exposed on the surface due to weathering and erosion (Goldman 1933;Bruthans et al. 2009). The overwhelming majority of reported salt diapirs in the Gulf of Mexico and elsewhere are found at considerable depths (> 160 m) and less frequently at surface or near-surface elevations (Escher and Kuenen 1928;Beckman et al. 1990). The very few salt diapirs that exist on the surface are found in hyperarid areas (e.g., Zagros Mountains: Alsop et al. 2016; Farasan Islands and Jazan city in Saudi Arabia: Almalki et al. 2015), where the arid conditions promote their preservation. This case study summarizes the deformation rates above one of the few well-exposed diapirs in the city of Jazan in southwest Saudi Arabia.
Over the years, salt diapir studies have received attention from the geologic community, given that most of the world's hydrocarbons reside in salt basins (e.g., Gulf of Mexico, Lower Congo Basin, Persian Gulf, North Sea Campos Basin, and Pricaspian Basin) (Hudec and Jackson 2011;Wu et al. 2015;Jackson and Hudec 2017). Moreover, salt diapirs are being investigated as potential places for underground storage of hazardous material (Koyi 2001;Jackson and Hudec 2017).
The Gulf of Mexico has one of the world's largest salt diapir occurrences, where hundreds of salt diapirs were discovered onshore and under the gulf itself (Beckman et al. 1990). There has been increasing interest in identifying the distribution of salt diapirs in the subsurface using seismic reflection and refraction surveys (Ladd et al. 1976;Black and Voigt 1982) and in studying their modes and mechanisms of emplacement given their growth rates, structure styles, and role in forming hydrocarbon traps for oil and gas (Ewing and Ewing 1962;May and Hron 1978;May and Covey 1983;Wu et al. 2015;Amin and Deriche 2015). Salt diapirs have also been investigated using gravity and magnetic surveys (Peters and Dugan 1945;Nettleton 1962;Blank et al. 1987; Barbosa and Silva 2011) where salt diapirs will usually, but not always, show a gravity minimum (Chapman 1983). In gravity and magnetic surveys, one or more post filters are applied to enhance the contrast in the geophysical signal between the subsurface target (the salt diapir, in our case) and its surroundings (Almalki et al. 2014;Fairhead et al. 2011;Miller and Singh 1994;Oruc 2010;Verduzco et al. 2004). Filtering methods have been used to separate the longwavelength anomalies (regional), which are associated with deep/extensive subsurface 1 3 features, from short-wavelength anomalies (residual), which are associated with shallow/ small features (Gupta and Ramani 1980;Mickus et al. 1991). Tilt derivative filters (TDR; Verduzco et al. 2004) were successfully applied to aeromagnetic and gravity data over the Farasan Islands to map the general distribution of near-surface diapiric structures (Almalki et al. 2014(Almalki et al. , 2015. Near-surface diapirs could also be identified in low-conductive environments using ground penetrating radar (GPR; Siever and Elsen 2010;Gundelach et al. 2012). Forward modeling of integrated geophysical data (e.g., gravity and seismic data) can illustrate the spatial distribution of diapiric structures (Almalki et al. 2015).
These methods, although effective in many environmental and geological settings, have their limitations. Seismic reflection/refraction surveys are challenging to implement and costly if deployed over large areas. Aeromagnetic and gravity surveys are subject to loss of target resolution given the relatively large station intervals (Issawy et al. 2011), and GPR signals are attenuated in clay-rich areas or those having saline groundwater (Annan 2005).
The advent of the synthetic aperture radar (SAR) system that operates in the microwave frequencies of the electromagnetic spectrum has enabled precise monitoring of the deformation of the Earth's surface (Rosen et al. 2000;Bürgmann et al. 2002). SAR gets its name from the repeated radar pulses as it moves along its flight path creating a synthesized antenna length, several kilometers in length (Rosen et al. 2000;Simons and Rosen 2015). Interferometric synthetic aperture radar (InSAR) refers to SAR applications involving combining signals from two SAR antennas. Ground motion can be measured by imaging the surface at multiple times usually months to years over the same location (Rosen et al. 2000). InSAR has two main components, phase and amplitude. Amplitude provides information on the intensity of the backscatter signal and the phase on the distance between the sensor and the target (Balmer and Hartl 2010;Simons and Rosen 2015). Two methods are commonly used to measure the phase, differential InSAR (DInSAR) and repeat-pass InSAR. The former has been used to monitor earthquake-and volcano-related deformation using before and after images, and the latter detects cm-to mm-scale deformation over extended periods of time (Bürgmann et al. 2002;Balmer and Hartl 2010).
Interferometric synthetic aperture radar (InSAR) provides cost-effective continuous measurements of land deformation over large areas with high precision (Gabriel et al. 1989;Devanthéry et al. 2016). InSAR has been used successfully to monitor subsidence associated with sediment compaction in deltas (Becker and Sultan 2009;Gebremichael et al. 2018;Delgado Blasco et al. 2019) or groundwater extraction (Chaussard et al. 2014;Othman et al. 2018), as well as uplift associated with volcanic activities (Bürgmann et al. 2002;Hooper et al. 2004Hooper et al. , 2007. However, relatively few studies have applied InSAR to monitor displacement and deformation associated with salt intrusions in arid or rural environments. These include the Precambrian Hormuz salt intrusions in the Zagros Mountains of southern Iran (Aftabi et al. 2010) and Tertiary salt intrusions in the Great Kavir basin in northern Iran (Baikpour et al. 2010). To the best of our knowledge, no published studies to date have applied InSAR technologies to assess the environmental impacts associated with salt intrusions in urban settings. In this study, we apply InSAR technologies to examine the deformation rates associated with the rise of the JZD in the Jazan city.
Permanent GPS stations are used to monitor crustal deformation and movements on a global scale (Abdrakhmatov et al. 1996;Becker 2015, 2019;Zerbini et al. 2017).
In this study, we compare the radar-based deformation rates to those extracted from GPS data (Global Positioning Systems) to evaluate and compare the radar-based velocities against the GPS-based absolute deformation rates. In conducting this exercise, we identified the location of the GPS station on the radar scenes and extracted the deformation rates over the GPS station from all datasets in order to directly compare the results (Simons 1 3 and Rosen 2015). This step enables the generation of a deformation rate product that takes advantage of the large aerial coverage of the radar-based scenes and the high accuracy of the GPS deformation rates.
Passive seismic data can detect the physical contrast between layers with varying shear wave velocities and densities and were thus used in this study to help identify the extent and depth of the salt layer over the JZD and its immediate surroundings. Microseismic tremors have been measured by seismologists and are generally regarded as background noise (Bonnefoy-Claudet et al. 2006). More recently, researchers have realized the significance of using the 'background noise' in geological investigations. The horizontal-to-vertical spectral ratio (HVSR) method was used to identify the depth to the bedrock underlying the unconsolidated layer (Delgado et al. 2000). To date, this passive seismic technique has been successful in examining a wide range of geophysical and geological problems including seismic hazards (Nakamura 1989) and mapping sediment thickness (Seht and Wohlenberg 1999;Chandler and Lively 2014;Johnson and Lane 2016).
The gravity method is categorized as a potential field geophysical method since it measures naturally existing magnetic and gravity fields created by the Earth. Recent advances are enabling precise measurements of the gradients of the gravity field (Lee 2001;Hammond and Murphy 2003) that could potentially be used in subsurface exploration. Gravity techniques identify areas of contrasting density in the subsurface by collecting surface measurements that detect variations in the Earth's gravitational field. For many conceivable sources (e.g., salt bodies, water-filled cavities, soft sediment-filled cavities, or partially filled cavities), the density contrast with the surrounding rocks shows as a low relative anomaly on the Bouguer map. Several gravity studies have been conducted in the areas surrounding the JZD. A regional gravity survey (Gettings 1977(Gettings , 1983) provided a low-resolution Bouguer anomaly map over the entire Jazan Province, and a high-resolution 50 m gravity survey was conducted over the Farasan Islands (Almalki et al. 2015). In this study, we report the first high-resolution 1 km gravity survey over the JZD and its surroundings.
We apply an integrated approach that takes advantage of InSAR technologies, gravity measurements, and passive seismic data to accomplish the following: (1) identify the distribution of salt diapirs at surface or near-surface elevations, (2) investigate the nature and rates of deformation associated with the intrusion of these shallow diapirs, and (3) examine their potential role as an environmental hazard. In this study, we employ the Jazan diapir (JZD) in the city of Jazan within the Jazan Province, Saudi Arabia, as a test site for our novel integrated approach. The JZD intrudes a thick (5 km) sequence of Tertiary and Quaternary deposits flooring the coastal plain (Blank et al. 1987). Ongoing ground deformation associated with salt diapir ascent has caused the destruction of infrastructure, and residential and governmental buildings (Erol and Dhowian 1988;Fatani and Khan 1993).

Geologic, Hydrologic, and Climatic Settings
The city of Jazan is a main port on the Red Sea in southwest Saudi Arabia. The city is located some 60 km north of the border of Saudi Arabia with Yemen ( Fig. 1a, b). Since the 1990s, the old city of Jazan has been witnessing urban expansion along the coastline to its north, south, and east. The city of Jazan occupied an area of 262 km 2 in 2008, 354 km 2 in 2014, and 409 km 2 in 2018, an increase in size amounting to 60% between 2008(Al-Zubieri et al. 2018Abd El-Hamid et al. 2019). A part of the coastal urban expansion was on reclaimed coastal regions and active sabkhas to the west of the old city.
The city lies within a wide (50 km) coastal plain that is bound by the Red Sea Hills to the east and by the Red Sea shoreline to the west (Fig. 1a). Valleys (wadis) crosscut the coastal plain and channel runoff from precipitation over the Red Sea Hills toward the coastal plain. These wadis are floored by thick alluvial and floodplain terrace deposits ). The coastline is decorated by eolian sand and sabkha deposits. The latter vary in thickness (10-16 m thick) and consist of sabkha crust (1-2 m), sabkha complex (5-9 m), and sabkha base (4-5 m) (Erol 1989;Dhowian 2017); the bottom two layers are characterized by fine clay deposits and silty sandstone (Erol 1989;Youssef et al. 2011;Youssef and Maerz 2013) (Fig. 1a). The coastline deposits and the coastal plain are underlain by thick (up to 5 km) Quaternary and Tertiary clastic deposits that overlie basement blocks that were down-dropped along northwest-trending extensional  Erol [1989] and Gillman [1968]) showing the extent of the study area, the distribution of lithological units, the JZD known extent, areas where land deformation was reported (light green polygons), borehole locations (BH1 and BH2), GPS and absolute gravity stations, and the locations where Tromino measurements were acquired. b The locations of the known salt diapirs including the JZD diapir in Jazan city, Samtah diapir to its south, Farasan Island diapirs to its west, and the Al Salif diapirs off the coast of Yemen. c Radar scene footprints showing all three orientations of SAR data used including Envisat descending track 6 (red polygon), Sentinel-1 ascending track 116 (green polygon), and Sentinel-1 descending track 6 (blue polygon) faults during the Red Sea opening (Blank et al. 1987). The northwest-trending faults crosscut the earlier northeast-trending Precambrian fault systems (Basahel et al. 1983).
Several diapirs have been reported within, or proximal to, the Red Sea coastal plain, one of which crops out within Jazan city, the JZD, which is the focus of this study (Fig. 1). The JZD covers an area of 2 km 2 , rises some 50 m above the flat coastal plain, and extends in a northwest-southeast direction (Gillmann 1968). Approximately 30 km to the southeast of the JZD is a northwest-southeast-trending 1 km 2 diapir proximal to the city of Samtah ( Fig. 1b; Blank et al. 1987). Some 100 km to the south in Yemen, a massive salt diapir (Al Salif diapirs) crops out over an area of 6 km 2 , extends for 4 km in a northwest-southeast direction, and rises some 45 m over the coastal plain ( Fig. 1b) (El-Anbaawy et al. 1992;Davison et al. 1994Davison et al. , 1996Bosence et al. 1998). Offshore of the JZD and some 50 km to the west are the diapirs of the Farasan Islands (Fig. 1b). These are circular to elliptical bodies that range from 3 to 35 km in diameter, extend in a northwest-southeast direction, and underlie an uplifted coral reef deposit (Almalki et al. 2014). Extensive field and high-resolution gravity data indicate that salt is present in some areas within the islands at depths as shallow as some 20 m (Almalki and Bantan 2016). The parent salt body that makes up the above-mentioned diapirs is believed to be of Miocene age (Blank et al. 1987;Bosence et al. 1998) and could have been mobilized by the onset of rifting along the Red Sea (Almalki et al. 2015).
The JZD is formed of halite with a caprock consisting of a chaotic mixture of steeply dipping gypsum and anhydrite of the Baid Formation (shale and sandstone) that makes up the overburden sediments (Schmidt and Hadley 1985). The succession, depth, and thickness of lithologic units that were intercepted in a borehole north of the JZD and which were intruded by the JZD diapir from top to bottom are: medium-to fine-grained sandstone (0-200 m); anhydrite and gypsum with argillaceous silt and organic-rich shale (200-400 m); graywacke and polygenic conglomerate (400-450 m); anhydrite and gypsum (450-850 m); silty sandstone, gray-green shale, and shale-rich organic matter (850-900 m), followed by a massive halite core (Gillmann 1968). Across the JZD, the overburden varies in thickness due to the differential weathering, uneven distribution of the caprock and localized dissolution of the halite layer. Several dissolution sinkholes are located within the known outcrop with diameters ranging from 4 to 6 m (Basyoni and Aref 2015).
The alluvial aquifer is the primary source of groundwater in the coastal plain where the watersheds draining the Red Sea Hills collect precipitation and channel the accumulated runoff toward the Red Sea coastal plain aquifers as surface runoff and/or groundwater flow. The precipitation over the coastal plain and the surface runoff infiltrate and recharge the shallow alluvial aquifer, as does the groundwater flow from the fractured basement aquifer (Sultan et al. 2019;Alshehri et al. 2020).
Humidity is relatively high (45-65%; Basyoni and Aref 2015), the average temperature is 23 °C (~ 73 °F), and rainfall is irregular and usually occurs over relatively short, heavy storms. Average annual precipitation (AAP) over the coastal plain is low (AAP: 100-500 mm) compared to that over the Red Sea Hills to the east (AAP: 700 mm; Batayneh et al. 2012;Youssef and Maerz 2013).

Data and Methods
Our integrated methodology encompassed the following steps: (1) compile and inspect reported damages for buildings and infrastructure on or surrounding the JZD; (2) conduct radar interferometric studies to assess the land deformation over the JZD and its surroundings; (3) extract the average vertical component velocity from the GPS data and compare these velocities to those obtained from the InSAR products; (4) conduct a gravity survey over the JZD to constrain the nature of its contacts with the surrounding country rocks; (5) collect passive seismic data over the JZD and its surroundings to determine the thickness of overburden overlying the salt diapir, and (6) conduct limited drilling to examine inferences from the above-mentioned geophysical methods.

Field Campaign Observations
Two field campaigns (2016 and 2018) were conducted to inspect the distribution of reported damages in buildings and infrastructures over the JZD outcrop and in its surroundings. Figure 1a provides the distribution of eight zones identified within the JZD outcrop in which building and infrastructural cracks, defects, damages and failures were reported in the 1980s and 2000s (Erol and Dhowian 1988;Fatani and Khan 1993;Youssef et al. 2011;Youssef and Maerz 2013). During the field campaigns, geophysical data (gravity and passive seismic) and borehole data were collected. Spatial correlations were also conducted between the distribution of reported and observed structural damages and the radar-based distribution of areas of high deformation (subsidence and/or uplift).  Fig. 1c were used to assess land deformation and identify variations in deformation patterns and intensity through time. Two methods were applied, persistent scatterers interferometry (PSI) and the small baseline subset (SBAS) (Figs. 2, 3). Both methods were used to test for consistency of the extracted land deformation velocities (Hu et al. 2019). Line-of-sight (LOS) velocities were extracted from both Sentinel-1 and Envisat datasets, but SBAS LOS velocities were extracted from Sentinel-1 data only due to the large spatial and temporal baselines in the available Envisat dataset. Two different software algorithms were used to process the above-mentioned radar datasets. The Envisat scenes were processed using the Stanford Method for Persistent Scatterers (StaMPS) algorithm (Hooper et al. 2004(Hooper et al. , 2007, while SARscape 5.5 (https ://www. sarma p.ch/wp/index .php/softw are/sarsc ape/) was used to process both of the Sentinel-1 datasets. The Sentinel-1 results were processed to investigate ongoing patterns of deformation and to compare the extracted velocities and patterns to those from the Envisat dataset for the earlier period (2003)(2004)(2005)(2006)(2007)(2008)(2009). Correlations between the extracted InSAR displacement using the SBAS and PSI techniques were examined over areas of interest and along selected profiles (Fig. 4).

Envisat PSI Processing
Seven level 1 SLC (single look complex) descending synthetic aperture radar (SAR) scenes (track 6) covering the period from 2003-2009 were acquired by Envisat's Advanced Synthetic Aperture Radar (ASAR) sensor. Although the number of available scenes was small (Fig. 2a), it was sufficient to estimate the deformation patterns in and around the diapir given that the study area is lacking a vegetative cover that could reduce the overall coherence of the individual pixels (Barrett 2012;Catry et al. 2016).
InSAR PSI was applied to a stack of Envisat scenes to quantify the deformation over the above-stated period. The PSI method restricts the phase unwrapping and analysis to coherent pixels, termed persistent scatterers (Ferretti et al. 2001). The StaMPS algorithm was used to extract/identify such pixels that exhibit phase stability in the study area (e.g., buildings, outcrops, rocks, dwellings, or utility poles). The StaMPS method was selected due to the algorithm's capability to retain low-amplitude pixels with poor signal-to-noise ratio but that exhibits phase stability over the investigated interval (Crosetto et al. 2015) despite poor temporal decorrelation between scenes. This is the reason behind the selection of the StaMPS technique over SARscape for the analysis of the limited number of available Envisat scenes.
The SAR images were co-registered and realigned to a selected primary scene (05/10/2004; Fig. 2a). The primary scene was selected based on minimizing overall spatial and temporal baselines. Two-pass differential interferograms were then generated between the primary, and each of the secondary reference scenes. Pixels that rendered phase stability were identified as persistent scatterers (PS). A coherence threshold of 0.3 was selected to identify PS points (Hooper et al. 2004). The topographic phase contribution was simulated using the 1 arcsecond (30 m) Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) and removed from the overall deformation estimate. Delft precise orbits were used to refine orbital information. Due to the limited number of acquisitions Fig. 2 Spatial baseline-time plots for Envisat and Sentinel-1 tracks. a descending Envisat track 6 PSI (7 scenes). b ascending Sentinel-1 track 116 PSI (83 scenes). c descending Sentinel-1 track 6 SBAS (32 scenes). The figure also shows the distribution of the primary scenes (yellow diamonds), the reference SAR images (green diamonds), and the interferometric pairs in blue with very long and irregular temporal baselines, the SBAS techniques were not applied to the Envisat stack.

Sentinel PSI Processing
Eighty-three ascending Interferometric Wide (IW) mode scenes (2014-2018) were acquired from the ESA C-band (5.56 cm wavelength) Sentinel-1 satellite (track 116). The scenes were downloaded from ASF Data Search Vertex. The Sentinel-1 mission is a currently active constellation of two twin C-band satellites (Sentinel-1A and Sentinel-1B) operated by ESA (Torres et al. 2012;Potin et al. 2016). Same-track Single Look Complex (SLC) Sentinel-1 frames were mosaicked to cover the area surrounding the diapir in Fig. 3 Line-of-sight (LOS) land deformation (in mm/year) for the JZD and surroundings extracted from the Envisat and Sentinel-1 datasets plotted over a DEM (TanDEM-X ~ 12.5 m resolution) from the German Aerospace Centre (DLR) TerraSAR-X/TanDEM-X Satellite. Lighter colors indicate higher elevations within the diapir and darker colors indicate low to below sea-level elevations for the surrounding sabkhas. Negative velocities (orange to red) represent downward ground motion (subsidence), whereas positive (green) velocities represent upward ground motion (uplift). For visualization purposes, velocities ranging from − 2 to 2 mm/year were assigned the same symbology (beige) to highlight the rapid rates of deformation. a LOS land deformation extracted from Envisat ASAR SLC descending (track 6) scenes that were acquired between 2003 and 2009 and processed using StaMPS. b LOS land deformation extracted from Sentinel-1 ascending (track 116) scenes that were acquired between 2014 and 2018 and were processed using SARscape. c LOS land deformation extracted Sentinel-1 SBAS descending (track 6) scenes that were acquired between 2014 and 2018 and processed using SARscape. Land deformation velocities along northsouth transects (a-a'; b-b'), east-west transects (c-c'; d-d') and northwest-southeast (e-e') transect are plotted in Fig. 4. Also shown is the location of the GPS station, the known extent of the JZD (white polygon) and five areas in the immediate northern, southern, and western surroundings of the known extent of JZD (outlined by red polygons) that are undergoing uplift and are susceptible to future land deformation the ascending geometry. The original scenes were cropped to a smaller area covering the diapir and a large portion of the surrounding coastline, approximately 50 km to the north and south of the JZD diapir outcrop. Both ascending and descending geometries were processed. The ascending geometry was selected for having less atmospheric noise and better coherence over the study area. ESA's Sentinel-1 precise orbits were used to refine the orbital information. Available ALOS World 3D DEM (30 m resolution) was used to remove the topographical phase component. Stacked scenes were processed using SARscape at a ground resolution of 60 m. A single primary scene was selected (07/11/2017; Fig. 2b), the atmospheric phase component was estimated and removed by low-pass spatial (10 km) and high-pass temporal filtering (365 days; Ferretti et al. 2000). A coherence threshold of 0.4 was selected for the final geocoding step (Fig. 3b).

Sentinel SBAS Processing
Thirty-two descending Sentinel-1 (track 6) IW mode scenes (2014-2018) were processed using the SBAS method. SBAS applications rely on (1) creating a stack of interferograms with small temporal and orbital separation (baseline) and (2) filtering out the topographic artifacts and atmospheric phase component using spatial and temporal information available in the processed data (Berardino et al. 2002). Sixty-three scenes were initially downloaded from the ASF Data Search Vertex. Thirty-one scenes were excluded from further processing due to their long temporal or spatial baselines or consistent poor interferograms Fig. 4 Distance (m) versus InSAR displacement velocities (mm/yr) plots derived from the PSI applications using Envisat (orange line) and Sentinel-1 (blue line), and SBAS Sentinel-1 (gray line) extracted along north-south ( Fig. 3: a-a'; b-b'), east-west (Fig. 3: c-c'; d-d'), and northwest-southeast (Fig. 3: e-e') transects. The known extent of the diapir is shown in dashed vertical lines with all of the remaining scenes in the stack. This could be related to one or more of the following: low signal-to-noise, low coherence, or atmospheric errors. The scenes were clipped to process areas immediately surrounding Jazan city. Stacked images were processed using SARscape with a ground resolution of 15 m. SLC image stacks were coregistered, and a primary image (04/17/2017; Fig. 2c) was selected. Differential interferograms were created for image pairs having small temporal and orbit baselines. Phase unwrapping was performed on the interferograms using Delaunay triangulation (Costantini and Rosen 1999). Twenty ground control points (GCPs) were used for the descending stack for orbital refinement by estimating and removing the residual phase ramps on the unwrapped interferograms. The GCPs were selected from the locations where no deformation was observed or expected. The first inversion of the SBAS method calculates a mean velocity field and possible topographic phase residuals using the unwrapped phases. Topographic phase residuals were subtracted from the wrapped interferograms. In the second inversion step, time series of displacement were estimated by the singular value decomposition (SVD) inversion approach using the refined wrapped interferograms (Berardino et al. 2002). Atmospheric phase component was estimated and removed by low-pass spatial (1200 m) and high-pass temporal filtering (365 days; Ferretti et al. 2000). Final velocity and time series of deformation along with the deformation velocity precision and topographical residuals were retrieved and geocoded in geographic coordinates (Fig. 3c).

GPS Data Studies
The SAR-based rates of surface deformation were compared to those extracted from the continuous GPS (cGPS) data from the continuously operating reference stations (CORS) of the Saudi Ministry of Municipal and Rural Affairs (MOMRA) Real Time Network (MRTN). The cGPS station (GZAN) is located in the southwest corner of the known extent of the JZD (lat: 16.880755°N, Long 42.542957°E; Fig. 1a). Three-component (east-west, north-south, and up-down) daily deformation/velocity data were retrieved from the permanent station. Only the vertical velocity component was used to compare the relative land deformation against the SAR-based estimates. The station covers roughly the same period as the Sentinel-1 data (2015)(2016)(2017)(2018). The daily GPS data were analyzed using the GNSS-Inferred Positioning System and Orbit Analysis Simulation Software (GIPSY-OASIS) developed by the Jet Propulsion Laboratory (JPL) at NASA (Wu et al. 1986). The output of the GIPSY-OASIS software was first converted from km/yr to mm/yr and then subtracted from the estimated mean velocity of daily values. A three-point average was applied to minimize seasonal effects, and a trend line was extracted for comparisons with the extracted radar results. The average vertical velocity component was extracted from the GPS over a period of approximately four years (Fig. 5) and was compared to the LOS velocities that were estimated from each of the three InSAR datasets.
The PSI and SBAS LOS velocities over the selected GPS station were obtained by averaging those of the PSI and SBAS scatterers within a circle (radius: 100 m) centered over the station. The inverse weighted distance (IWD) method was applied to assign more weight to the points closer to the GPS station. We estimated the error associated with the calculated trend values of the PSI, SBAS, and GPS time series by using a statistical approach that accounts for trend errors in the time series (Scanlon et al. 2016). Monte Carlo simulations were performed by fitting trends and other terms for many (n = 10,000) synthetic datasets. The standard deviation of the extracted synthetic trends was interpreted as the trend error (Table 1).

Gravity Data
Standard procedures were followed in the collection of field gravity data and in applying data reduction (Reynolds 2011). We conducted a detailed gravity survey by collecting 345 gravity stations over the study area and surroundings using two Scintrex gravity meters. The CG-5 Autograv system is a microprocessor-based automated gravimeter that has a measurement range of over 7000 mGal and a gravity resolution of 0.001 mGal. The CG-5 system is based on a quartz spring, which reduces the instrumental drift error and is highly resistant to thermal changes. The reading time was set to 120 s, and the seismic filter was enabled to reject seismic noises. The station interval ranged from 1 to 2 km, with larger distances over inaccessible areas such as the areas occupied by military bases and airports along the coastline. Previous regional studies over the study area and its surroundings used spacings ranging from 5 to 20 km. Wide-spaced gravity stations in regular or semiregular grid may reduce the survey resolution and could potentially overlook the presence of small near-surface diapirs. Stations were chosen along 24 east-west trending profiles over the salt  diapir and its surroundings. Different locations were tested for selecting the most stable and accurate gravity base station for the survey. Repeated readings at the selected base station were collected at the beginning and end of the day to determine the instrumental drift value. The relative gravity dataset was tied to the QIZAN (16.90933°N, 42.54883°E) absolute base station with an observed gravity of 978,488.379 mGal. The gravity data were corrected for temporal variations (drift and tides) using the Geosoft Oasis montaj version 6.4 software. The coordinates of stations and their elevations were obtained by real time kinematics (RTK) GPS surveys, with an error of ± 15 mm in elevation and approximately ± 0.004 mGal of error in gravity values. The Free-Air correction was applied, and a density of 2.67 g/cm 3 (Murata 1993) was used to generate the Bouguer anomaly map of the study area using Oasis montaj software (Fig. 6).

Passive Seismic Data
Passive seismic data were collected over a period from 10 to 14 minutes for 17 stations in and around the known extent of the diapir (Figs. 1a, 7a). The location of each station was established using a handheld GPS unit. Passive seismic data were collected using a TZ3 Tromino, a three-component, broadband seismometer. The Tromino is a small, portable instrument with an operating range of 0.1-64 Hz (Okada and Suto 2003). It measures background motion (e.g., industrial machinery, road, rail traffic, earthquakes, wind, ocean waves, etc.) to determine the thickness of the overburden overlying the bedrock, in our case the overburden over the salt diapir. The three components (east-west, north-south, and vertical) were processed (in 20 s processing windows) using Micromed Grilla software; processing involved Fourier transformation, spectral analysis, combining the two horizontal spectra (east-west and north-south), dividing the horizontal spectrum by the vertical spectrum, averaging all the processing windows, and smoothing the H/V curve (Konno and Ohmachi 1998;Johnson and Lane 2016). Data containing noticeable noise bursts were removed. The remaining data were reprocessed to produce a smoother curve over the spectral range, and a sharper peak on the HVSR plot. The frequency and amplitude of the dominant peak were identified on the HVSR plot for each station. The frequency of the H/V peak can provide an estimate of overburden thickness. The thickness of the soft layer (overburden) was estimated from the expression: H = V s /(4×f 0 ), where H is the layer thickness, f 0 is the measured resonant frequency from the H/V plot, and V s is the shear wave velocity. Passive seismic data and known depths to consolidated layers from two boreholes BH1 ( Fig. 1a; lat: 16.877618°N, long: 42.546327°E) and BH2 ( Fig. 1a; lat: 16.887648°N, long: 42.563869°E) were used to calculate the shear wave velocities over and to the east of the diapir, respectively. A shear wave velocity of 140 m/s was calculated at BH1 (station over the diapir) where a salt layer was intercepted at a depth of 4 m and an f 0 value of 8.75 was extracted from the passive seismic data collected from the nearest station to the borehole. A velocity of 122 m/s was calculated over BH2 (sabkha east of diapir) where a consolidated layer was intercepted at 18 m and an f 0 value of 1.69 was extracted from the nearest passive seismic station. Our calculated shear wave velocities are similar to those estimated using multichannel analysis of surface waves over the JZD and surrounding sabkhas (Alhumimidi 2020). The estimated shear wave velocities were then used to estimate depth to the caprock or halite body for the stations over the diapir only (Table 2).

Fig. 6
a Bouguer gravity anomaly map generated from 345 gravity station measurements in and around Jazan city that were tied to the absolute gravity base station (QIZAN; 978,488.379 mGal) north of the JZD (red circle). The contours are in mGals with an interval spacing of 0.5 mGal. The lower gravity values over the diapir and its surroundings are consistent with the presence of a salt body, and the closely spaced contours on the eastern side are consistent with a steep contact of the diapir with its surroundings. Also shown are the locations of boreholes BH1 and BH2, the known extent of the JZD (dashed green and black polygon), and cross-section A-A' (Fig. 1a). The inset map shows the distribution of the gravity stations Fig. 7 a WorldView-3 true color composite image over the Jazan city showing the distribution of 17 passive seismic survey stations, borehole (BH1 and BH2) locations, the known extent of the JZD (black polygon) and cross section A-A' (Fig. 1a). Red dots denote locations of passive seismic stations within the known extent of the JZD and green dots for stations outside (south, east, and west of the JZD). b amplitude (H/V) versus frequency (Hz) spectral plots displaying sharp peak H/V amplitudes and higher frequencies (red lines) over the stations within the known extent of the diapir. c lower peak H/V amplitudes and lower frequencies (green lines) outside the known extent of the diapir over the sabkhas. The higher the H/V ratio and the higher the frequency, the shallower the contact with the bedrock (caprock or salt layer) and vice versa for the deeper contact. Station T6 ( Fig. 7c; dashed green line) is outside of the known extent of the diapir, but it has an H/V amplitude and frequency similar to those stations over the diapir (7b). Station T2 is near BH1 and station T58 is near BH2 1 3

Drilling
Borehole BH1 was drilled as part of the required procedures and regulations for laying down foundations for new constructions. A halite body was reported at a depth of 4 m. We drilled a borehole (BH2) some 500 m to the east of the known extent of the JZD outcrop. The borehole was drilled within the sabkhas surrounding the JZD and reached a depth of 60 m. The drilling was conducted using a B53 mud-rotary hydraulic mobile drill and a Tricone bit. Seventeen samples were collected at irregular-spaced intervals (4-10 m) to identify the lithological sequence of the penetrated layers. BH1 and BH2 were drilled close to the Tromino stations T2 and T58, respectively (Fig. 7a). Figure 5 shows the extracted vertical velocities of the permanent GPS station (GZAN) during the period (2015-2018) and an average vertical displacement of − 0.76 ± 0.64 mm/yr. We compared the vertical velocities from the GZAN station to the PSI and SBAS LOS velocities. Comparisons of the deformation rates from the Jazan GPS station (Fig. 5) and those in Fig. 3 show a good general correspondence between the GPS velocities ( − 0.76 ± 0.64 mm/yr) and those from the PSI (Envisat: 0.81 ± 0.08; and Sentinel-1: − 0.27 ± 0.19 mm/yr) and SBAS (0.17 ± 0.37 mm/yr) over and proximal to the GPS station (Table 1). The radar-based velocities are similar to the GPS velocities when the estimated uncertainties in both the GPS and radar-based velocities are taken into consideration (Table 1). Thus, no calibration for the radar-based velocities was conducted. Figures 1, 8 provide the distribution of eight zones within the known boundaries of the JZD outcrop in which building and infrastructure damages were reported in the 1980s and 2000s (Erol and Dhowian 1988;Fatani and Khan 1993;Youssef et al. 2011;Youssef and Maerz 2013). These reported defects were inspected and verified during our field trips in years 2016 and 2018. Figure 8 shows two general patterns of deformation, the JZD outcrops witnessing uplift, whereas the sabkhas to the east of the JZD are undergoing subsidence. The center of the JZD is outlined by a dashed white and black polygon; the areas subtended by the latter and the dashed red and black polygon that outlines the known extension of the diapir are hereafter referred to as the periphery of the diapir. The uplift rates (from the Sentinel-1 PSI results) within the center of the JZD outcrop exceed the uplift rates along its peripheries. The center of the diapir includes zones 3, 4, 7, and 8 and its periphery encompass zones 1, 2, 5, and 6. In general, the uplift rates for the center of the diapir range from 2 to 4 mm/yr and those for its periphery range from 1 to 2 mm/yr. One should not expect to observe these general patterns throughout the entire diapir. Deviations could result from the accumulation of precipitation in topographic lows which could lead to enhanced local dissolution of the salt layer. Also, there could be variations in the salt depths under the diapir from a few to tens of meters leading to disparities in land deformation rates over the selected zones . Figure 3 shows the spatial distribution of land deformation from three radar interferometric datasets, Envisat PSI Fig. 3a), Sentinel-1 PSI (2014Fig. 3b), andSentinel-1 SBAS (2014-2018;Fig. 3c). Inspection of Fig. 3 shows a general correspondence in the patterns and rates of deformation between the three products. For example, over a period of 15 years, the average uplift rates over the center of the diapir are similar in all products (PSI Envisat: 2.7 mm/yr; Sentinel-1: 2.1 mm/yr; SBAS: 1.8 mm/ yr), so are the average rates over the diapir periphery (PSI Envisat: 1.8 mm/yr; Sentinel-1: 1.3 mm/yr; SBAS: 1.7 mm/yr) and the sabkhas (PSI Envisat: − 1.7 mm/yr; Sentinel-1: − 1 mm/yr; SBAS: − 1.9 mm/yr). The sabkha areas are subtended by the dashed red and black polygon and the blue and white line immediately to the east of the JZD in Fig. 8. There are deviations from the above-mentioned overall patterns. Five areas ranging from 0.2 to 0.6 km 2 outside the known boundaries of the JZD (Fig. 3b) are undergoing uplift (average deformation rates: PSI Envisat: 1.8 mm/yr; Sentinel-1: 1.3 mm/yr; SBAS: 1.7 mm/yr) and have been identified in this study as future areas of concern.

Findings
The above-mentioned correspondence between the three datasets could be also observed along north-south (a-a'; b-b'), east-west (c-c'; d-d'), and northwest-southeast (e-e') transects (Figs. 3, 4). The transects show compatibility between the radar interferometry datasets, both temporally and spatially. In other words, the patterns and rates of deformation have been fairly consistent over the Envisat and Sentinel-1 investigation periods. The transects reveal that the uplift (up to 4 mm/yr) is not only restricted to the  (Erol et al. 1988;Fatani et al. 1993) and verified during our field campaigns. The average deformation rates were calculated over the center, periphery, and the sabkha east of the diapir. The area identified as the center of the diapir is outlined by a dashed black and white polygon, and the area subtended by the latter and a dashed red and black polygon is the periphery. The areas east of the diapir that are subtended by the dashed red and black polygon and the blue and white line are occupied by sabkhas. Deformation rates extracted from Sentinel-1 PSI results are shown as contours in the background with a contour interval of 1 mm/yr. Green contours indicate uplift, and orange to red contours are areas witnessing subsidence. The background is a DEM (TanDEM-X ~ 12.5 m resolution) diapir boundaries but also extends to its immediate surroundings in the north, west, and south. This phenomenon is not observed on the eastern side. In contrast, subsidence (down to − 4 mm/yr) is prominent to the east of the diapir (Figs. 3, 4).
One potential explanation for the observed patterns of deformation is that the contact between the diapir and the surrounding rocks it intrudes is near-vertical on the eastern side, but less so on the northern, western, and southern boundaries. If we were to drill in the immediate surroundings of the diapir boundaries, we predict the diapir will be intercepted at shallow depths north, south, and west of the diapir, but not to its east. The significance of this finding, if true, is that continuous uplift of the diapir will in future likely compromise the integrity of the buildings and infrastructures to the north, south, and west of the diapir, but not to the east. Development to the east of the diapir should therefore suffer less damage over time and incur fewer ongoing remedial expenses due to diapir uplift.
Additional evidence for the postulated extension of the salt diapir in the subsurface and its boundaries with its surroundings comes from the terrestrial geophysical datasets (gravity and passive seismic) and drilling. The Bouguer anomaly map of the study area shows a low-gravity anomaly over the JZD and a rapid change in the gravity field from the west to the east direction which is represented by a narrow separation between contour lines (contour lines are closer to each other) and a more gradual change in the north, south, and west directions (Fig. 6). The low-gravity anomaly over the JZD reflects the density contrast between the JZD salt body (~ 2.2 g/cm 3 ) and its surroundings (2.67 g/cm 3 ; Jackson and Talbot 1986;Warsitzka et al. 2013). The Bouguer anomaly over the eastern side of the salt diapir reveals a steep gradient that could be indicative of a steep boundary for the salt diapir, whereas the remaining diapir boundaries show a less steep dip with contour lines further apart. The lowest gravity values (7.5 mGal) are observed over the center of the diapir, and the diapir known boundaries are approximately marked by the 8.5 mGal contour. A rapid increase in Bouguer gravity anomaly values are observed to the east of the diapir (up to 12 mGal), but not to the north and south (up to 9.5 mGal).
Passive seismic data were collected inside and outside the known boundaries of the JZD (Fig. 7a). The peaks in the HVSR are caused by shear wave impedance contrasts between overburden sediments and the caprock, salt layer, or other compact lithologic units that could be interlayered with the soft sediments. The HVSR peaks over, and to the west of, the diapir outcrop yield a high frequency (5-10 Hz) and a strong and clear HVSR peak indicative of shallow depth to the salt or caprock layer (Fig. 7b). To the immediate east, the peaks are not as strong and the frequency is low (1.5-3 Hz), possibly indicating a flank of the diapir with a steep dip and the presence of a compact sand layer (not salt or caprock). The low-amplitude resonance peak may indicate a weaker acoustic impedance contrast across the boundary (Fig. 7c). Two boreholes have been drilled to verify the presence of a shallow salt body consistent with Tromino stations T2 (BH1, over the diapir) and T58 (BH2, outside the diapir outcrop to the immediate east) (Figs. 1a, 7a). Only BH1 with a depth of ~ 4 m to the halite body was used to calibrate the Tromino stations over the diapir (Fig. 7). Depths calculated for the other Tromino stations can be found in Table 2. The one exception was Tromino station T6, which showed a similar HVSR plot and peak to that of the stations over the diapir. When calculated, the depth at T6 is ~ 6.6 m, very similar to those calculated in Table 2. Thus, this station is most likely over a shallow halite body or caprock. The borehole near T58, BH2 was drilled to ~ 60 m, and no salt was intercepted to that depth, and therefore no stations outside of the diapir were calibrated.

Conclusions
An integrated geophysical and remote sensing-based application was applied to delineate and quantify land deformation associated with the emplacement of the Miocene JZD within the Jazan city, Jazan Province, Saudi Arabia. Radar interferometric techniques (PSI and SBAS) were used to identify areas undergoing land deformation and to quantify the nature (uplift versus subsidence) and rates of deformation over the JZD and its surroundings. Two methods and three InSAR datasets (PSI: Envisat and Sentinel-1; SBAS: Sentinel-1) acquired over the JZD and its surroundings covering a period of 15 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) were utilized to monitor spatial and temporal land deformation over the diapir and its surroundings, examine the consistency of results, and identify changes in patterns and rates of deformation through time.
Analysis of the three InSAR datasets reveals the following: (1) similar spatial patterns and rates of deformation, where the center and the periphery of the diapir are witnessing uplift at rates of up to 4.7 mm/yr and 4.4 mm/yr, respectively, whereas the sabkhas to the east are subsiding at rates of up to − 7.5 mm/yr; (2) the ongoing deformation patterns persisted throughout the periods covered by the investigated Envisat (2003)(2004)(2005)(2006)(2007)(2008)(2009) and Sentinel-1 (2014Sentinel-1 ( -2018 operational periods, suggesting that the deformation will most likely continue in the future; (3) the areas where most of the deformation and damage of buildings and infrastructure have occurred correlate spatially with areas of high uplift rate within, or proximal to, the center of the diapir; and (4) five areas in the immediate northern, western, and southern surroundings of the JZD are undergoing uplift (up to 3.7 mm/ yr), in contrast to the subsiding sabkhas to the east. These observations were interpreted to indicate that the salt diapir in general and the areas in its immediate surroundings, with the exception of the eastern side, have been undergoing uplift throughout the past two decades and are projected to continue in the near future.
One way to explain the absence of uplift on the eastern side of the diapir is the presence of a steep contact on the eastern side of the JZD compared to its other sides. This suggestion is supported by gravity data that show a steep gradient probably indicative of a steep boundary for the salt diapir, whereas the remaining diapir boundaries show a gentle slope with contour lines further apart. This suggestion is also consistent with the findings from passive seismic and drilling data. The HVSR peaks over, and to the west of, the diapir outcrop yield a high frequency (5-10 Hz) and a strong and clear HVSR peak, indicative of shallow depth to the salt or caprock layer. To the immediate east, the peaks are not as strong and the frequency is low (1.5-3 Hz), possibly indicating a flank of the diapir with a steep dip. A salt layer was intercepted at a depth of 4 m over the western margin of the diapir but none was intercepted in a borehole east of the diapir to a depth of 60 m (Fig. 1a).
We suggest that additional near-surface diapirs in the vicinity of the JZD could potentially be identified using the geophysical and remote sensing-based techniques applied in this study. The application of these methods would help in assessing the potential threat of diapir rise in the region and would greatly assist in informing local planning and development efforts. Moreover, such investigations could direct attention to possible areas of interest or concern that may be threatened by ongoing deformation and help preserve life and property in regions similarly threatened elsewhere in the world.