Nationwide Subsidence in Bahrain Island: Drivers and Implications for Relative Sea-Level Rise

The IPCC-AR6 report suggests that the sea-level in the Persian Gulf is expected to rise by 2.1–4.9 mm/year by 2100, which is lower than the global projections of 4–14 mm/year. However, a central nationwide ground subsidence in Bahrain's low-lying island can aggravate these figures. The island is only 14 km wide, yet we observe a 7 × 12-km ground subsidence with vertical rates up to 20 mm/year for 2016–2021. We forecast a subsidence of up to 1.5 m by 2100 by extrapolating the average rate. To characterize the subsidence, we use the following methods: i) SBAS-InSAR with 94 and 131 C-band Sentinel-1A orbital radar acquisitions in ascending and descending directions for 2016–2021, respectively, and ii) Stacking-InSAR with eight L-band ALOS-1-PALSAR acquisitions for 2007–2010 and six ALOS-2-PALSAR-ScanSAR acquisitions for 2016–2021. We suggest three causes associated with this subsidence: First, exploitation of the Awali oil field; second, local aquifer depletion; and third, subsurface dissolution of anhydrites and chalky limestones. To assess the potential impact of subsidence on the shoreline, we measure the coastline evolution at three undisturbed beaches on the West coast from 1985 to 2021, excluding areas that underwent land reclamation. Using subpixel shoreline detection analysis from 308 Sentinel-2 and Landsat (L5, L7, and L8) acquisitions, we observe that the selected shores remain stable. However, observations reveal a shoreline retreat of up to 5 m/year on the southwestern coast for 2003–2014, surpassing the estimated rate of 0.08 m/year derived from tide gauge data. This nationwide subsidence should be considered when forecasting coastal infrastructure planning in Bahrain. We recommend performing a similar analysis in other low-lying Gulf islands, where oil exploitation occurs.


I. INTRODUCTION
G ROUND deformation mapping derived from radar inter- ferometry (InSAR) started in the 1990s [1], and today, the amount of SAR images available allows to produce ground displacement time series using algorithms, such as the Small BAseline Subset (SBAS) [2].Ground deformation can have various origins, including compaction of sedimentary rocks hundreds of meters deep in petroleum reservoirs [3].In that case, monitoring ground surface deformation over oil fields, along with semi-infinite elastic modeling, helps understanding the dynamic of petroleum reservoirs and supports complex withdrawal activities [3], [4].Moreover, collapses at the surface and subsurface can result from active karstification [5] and mining [6].In addition, land subsidence can also occur in overexploited confined aquifers, where a decrease in pore pressure and an increase in effective stress [7] lead to the rearrangement of compressible sediments, such as silts and clays.Similar deformations are observed with InSAR in populated alluvial valleys, such as several valleys in Mexico [8], Las Vegas valley [9], or the central valley in California [10], but also in populated deltaic plains, such as Jakarta [11], [12].When subsidence occurs on the shore, the sea-level increases in relation to the land.This phenomenon, known as relative sea-level rise (SLR) [13], is widely observed with InSAR, e.g., in the Mekong Delta, Vietnam [14], in California [15], on the northern coast of the Gulf of Mexico [16], in the Korean peninsula [17], and in Tianjin, northern China [18].
Bahrain's main island, one of over 100 small low-lying islands (<15 km wide) distributed throughout the Persian Gulf [19], is a notable example of an island hosting a historic oil field.Other islands of the Gulf, such as Kish island [20] and Qaruh island [21], are also located over or near hydrocarbon reservoirs.In Bahrain, the crude oil production (including the natural gas liquids and additives extracted from the ground [22]) increased by 47% between mid-2010 and mid-2011 (see Fig. 1).Considering the significance of the extracted fluid volumes from the subsurface, we hypothesize that it can generate measurable ground deformation, similar to other studies, such as the land subsidence observed over the Belridge and Lost Hills oil fields in California [23].As such, this accelerated fluid extraction between 2010 and 2011, can be plausibly associated with a sharp increase in subsidence.These high extraction rates, associated with the increase in energy demand, initiated our interest in assessing the potential occurrence of ground deformations over the main oil field in Bahrain to evaluate its implications on relative SLR and shoreline retreat as a case study representative of several oil and gas fields in islands in the Gulf.It is important to note that our SAR dataset has a gap from 2010 to 2016 to cover the beginning of this potential accelerated subsidence.
Climate change has adverse impacts on the sustainability of small islands developing states (SIDS) [24], such as Bahrain's islands [25], [26].For instance, the effects of the SLR constitute a significant submersion threat to SIDS and low-elevation coastal zones (LECZ) [27].Moreover, groundwater freshness and low-lying arable soil are vulnerable to seawater intrusion arising from SLR [28], [29].Furthermore, the increase in storm surges damages coastal infrastructures [30].
The sea-level in the Persian Gulf varies as a function of regional [31] and local climate variabilities.In the gulf, the tide exhibits a semiannual component associated with local meteorological conditions; therefore the corresponding sea-level variation produces higher levels in summer than in winter.For instance, in 2000, the meteorological station in Dhahran, Saudi Arabia, located ∼20 km away from Bahrain, recorded a 23 cm increase in sea-level during the summer season compared to winter [32].High tides combined with low atmospheric pressure and strong winds can create surges [33], which could lead to coastal flooding when they occur episodically over vulnerable LECZ.Moreover, Cabanes et al. [34] suggest that hydrographic conditions, particularly water thermal expansion, also impact the sealevel variation leading to an overestimation of the 20th-century global SLR (GSLR) measured at tide gauges.However, Miller and Douglas [35] refute the authors' conclusions suggesting that the effect of temperature and salinity is negligible compared with the influence of mass changes, which is a primary driver.
The tide gauges network used to monitor the SLR must be representative of the true sea-level and should not be subject to any factors that could interfere with the measurements, such as vertical land motion (VLM), i.e., several tide gauges sink along with the shore in the Persian Gulf [38].When VLM occurs, the observed SLR should be corrected.For Bahrain, we combine the tide gauge rate located in the city of Manama [38] with the corresponding VLM [38], [39] (see details in Appendix A) to obtain a corrected-absolute SLR rate ranging from 2.77(±0.56) to 2.95(±0.55)mm/year.Assuming a steady sediment transport flow (i.e., absence of beach accretion or erosion), this corrected rate would result in a cumulative SLR of 21.3-22.7 cm by the year 2100.Moreover, this absolute SLR would result in an average shoreline retreat of 85 mm/year in present-day conditions, and 60-140 mm/year after 2100 under the RCPs scenari (see details in Appendix B).Al-Jeneid et al. [29] also modeled the shoreline retreat in Bahrain, resulting from the predicted SLR in 2100, and mapped future potential coastal flooding areas; Hassan and Hassaan [40] performed the same analysis for the coast of Kuwait (North-West of the Persian Gulf).The government of Bahrain implemented strong measures to address these threats and increased its landmass.Thus, dredging and land reclamation policies were initiated in 1982 [41], followed by new guidelines in 2008 [42].In 2015, 110 km 2 or 14% of the total land of the kingdom of Bahrain was reclaimed from the sea [43].However, the consequent dredging increases the coastal turbidity, causing damage to the marine ecosystems, including seagrass beds, mangroves, coral reefs, and muddy shores [43].
Herein, we investigate the nationwide spanning subsidence in Bahrain, which has significant implications on the relative SLR, shoreline retreat, and overall nation's sustainability.Furthermore, understanding the rates and drivers of this subsidence provides a unique insight for forecasting the submersion of small low-lying islands in the Persian Gulf under relative SLR.The latter is crucial for understanding the potential upcoming changes in the exclusive economic zone (EEZ) constituting maritime borders [44], affecting the sovereignty over oil and gas fields.Bahrain is a representative case study for all of the above.

A. Geology and Oil Production
A salt anticline in the N-S direction (see Fig. 2) hosts the Awali onshore oilfield, which has been exploited since 1933 [46].The base of the salt swell is ∼5.5 km deep, and the fold is 48 km long and 18 km wide, dipping 2°on the East flank and 4°o n the West flank [46], [47].The crest of the episodic growing anticline is eroded and flattened.Fig. 3 shows the structural Fig. 3. Sequence stratigraphy.We adapted the top cross section from [54] and the bottom one from [49].We outline these two profiles in Fig. 2.
The dome shape of the anticline forms structural traps where hydrocarbons accumulate under folded impermeable rocks, such as shale or anhydrite [51], referred to as cap rock.Bugti et al. [51] suggest that fault traps could also occur in the Awali field.For instance, numerous faults occur in the Wasia group formations, which contain the Mauddud reservoir, i.e., the largest oil reservoir of the Awali field [52].The Awali field comprises 16 oil and 6 gas reservoirs of fluids varying from tarry oil in the Aruma zone to dry gas in the Khuff zone [46].These reservoirs overlap in sedimentary traps, i.e., carbonate, for the most but also sandstone [46].
The fluid dynamics in these reservoirs are complex.For instance, in the mature Mauddud reservoir, there are occurrences of vertical and horizontal oil migrations, and aquifer influx on the reservoir's edge, in response to gas injection [53].

B. Climate and Hydrogeology
The climate in Bahrain is hyperarid, with rare rainfalls cumulating precipitation of <80 mm/year, and the total natural discharge from the aquifers exceeds the recharge [54].Surface water being inexistent, Bahrain's water management relies on aquifers and desalination [54].In 2011, the total water need of 387 Mm 3 /year was met with desalination (54.1%), groundwater abstraction (36.5% or 141.3 Mm 3 /year), and treated wastewater (9.4%) [55].
Over the first few hundred meters, the anticlinal structure shows a sequence of Neogene, Dammam, Rus, and Umm Er Radhuma (UER) Formations.The Dammam Formation shows Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
a succession of limestones/dolomites and shales/marls.The Rus Formation comprises a fractured chalky dolomitic limestone intercalated with anhydrites and shales in its upper part [54].Most of these anhydrite had dissolved in the central and eastern parts of the island, followed by a collapse of the overlaying formations, which also affected the confining layer and started an upward groundwater migration from the Rus-UER to the Dammam aquifer [54].The lithology of the UER Formation is similar to the Rus Formation, which is composed of karstified dolomitic limestones showing fractures and high heterogeneous porosity [49].The first part of the Dammam aquifer is confined at the top by a 10-60-m-thick claystone from the Neogene Formation.Underneath, the second part is semiconfined by the 9-15-m-thick orange marl member of the Dammam Formation.This regional aquifer crops out in central Saudi Arabia, where the recharge occurs and discharges in the eastern Arabian Peninsula [54].In our region of study, below this Dammam aquifer, the Rus-UER aquifer is confined by the 8-20-m-thick shale beds and anhydrite and shale deposits [54].In the center of the island, where the Rus formation outcrops from erosion, the Rus-UER aquifer is under unconfined condition.The Dammam aquifer is also unconfined where the Khobar member outcrops due to erosion [50].
Fig. 3 shows the aquifer system that holds the Dammam and the Rus-UER aquifers.The Dammam aquifer has been overexploited since the 1970s, with a groundwater abstraction exceeding the safe yield of 110 Mm 3 , reaching a maximum of 255 Mm 3 in 1998 [55], [56].Because of its overexploitation, the head levels have declined, and the hydraulic gradient reversed [56].The consequent encroachment of adjacent seawater and the underlying brackish aquifer increased salinity in the aquifer [49].In 2002, the Dammam aquifer was storing freshwater, with total dissolved solids (TDS) ranging between 3500 and 10 000 mg/L depending on the location [54].Underneath, the Rus-UER aquifer holds a brackish water lens with a TDS ranging between 8000 and 15 000 mg/L in 2003 [49].
Besides the issues related to water withdrawal and management, aquifers also play a role in gas recovery.In the case of highly permeable formation, water may encroach into the gas zone, trap gas pockets at high pressure and then reduce the final production by 30%-85% [57].Thus, understanding the dynamics and interactions of oil/gas fluids with aquifers is vital for optimized resource extraction and forecasting its impact on the subsurface.

C. Datasets
Our investigation uses three datasets to assess the central subsidence in the Awali oil field.Using the first, we produce four LOS deformation maps from SAR interferometry (InSAR): two with the Sentinel-1A SAR satellite (ascending and descending orbits), one with ALOS-1 PALSAR (ascending orbit), and one with ALOS-2 PALSAR ScanSAR (ascending orbit) (see Table I).The second is a temporal subset of the first dataset to extract the vertical and horizontal components of the LOS deformation maps within their common period (see Table II).Using the third collection, we detect the shoreline variations using the COASTSAT toolkit (see Table III).Furthermore, we use the NASADEM digital elevation model (DEM) [58] (∼10 m vertical accuracy [59]) to correct the topographic phase component in the interferograms.

III. METHODS
Our investigation implements two InSAR geodetic methods, SBAS-InSAR and Stack-InSAR, and the COASTSAT toolkit to assess shoreline retreats.These three methods are well established and widely used; therefore, we briefly summarize each to support an understanding of the reported results and their limitations.
1) SBAS-InSAR: Using 94 C-band Sentinel-1A scenes acquired in the ascending orbit direction between 2016-06-03 and 2019-12-21, we produce 940 interferograms with a spatial resolution of 140 m (see Table I).To improve the reliability of our results, we mask areas of low interferometric coherence.Then, we integrate selected interferograms (see Table II) with the intermittent SBAS algorithm [2], [67] to retrieve LOS time-series displacement velocity and its acceleration, as expressed in (1), while also accounting for intermittent coherent targets.We apply the same SBAS-InSAR methodology to 131 Sentinel-1A scenes acquired along a descending orbit direction during 2017-03-21-2022-01-30 where D is the displacement at time t, K 1 is the linear term, in mm/year, corresponding to the displacement velocity; K 2 is the quadratic term, in mm/year 2 , corresponding to the displacement acceleration; and K 3 is the cubic term, in mm/year 3 , corresponding to the change in displacement acceleration.We opt for the above cubic SBAS deformation model, in (1), to fit the nonlinear observed InSAR displacements associated with potential accelerations arising from changes in subsidence rates, occurring over the extended observation period of our study (2016-2022).
InSAR displacements are expressed along the LOS direction.It is the accurate way of presenting the results without making assumptions about horizontal deformation.However, these LOS results are not comparable between orbital sensors emitting at different look angles.Assuming that the horizontal deformation is null, the vertical contribution of the LOS displacement is equal to the latter divided by the cosine of the radar incident angle [68].For instance, a LOS displacement of 1 cm becomes a vertical displacement of 1.39 cm when the incidence angle (θ) is 44.0°(Sentinel-1A) and 1.21 cm when the angle is 34.3°(ALOS-1), which is a substantial difference of 18% between sensors for an identical vertical deformation.Such estimation is inaccurate when the horizontal deformation is not null because we do not consider the geometric error that varies as a function of the horizontal displacement and θ [69, p. 8-9].To remediate this vertical displacement inaccuracy, we combine two Sentinel-1A LOS velocities acquired with different geometry: one in ascending and the other in the descending direction, allowing us to extract their vertical and horizontal (eastward) contributions using a trigonometric relation [69].
2) Stacking-InSAR: From eight L-band ALOS-1 scenes acquired in the ascending direction between 2007-07-05 and 2010-07-13, we produce eight interferograms for the frames 500 and 510.We also produce seven interferograms from six L-band ALOS-2 ScanSAR scenes acquired in the ascending direction between 2016-08-08 and 2021-09-13 (see Table I).We mask areas of coherence below 0.3.Then, we remove the topographic phase component and correct the orbital inaccuracies for each interferogram.This number of interferograms is not large enough to fit a robust deformation model, as implemented in SBAS-InSAR.Instead, we perform the Stacking-InSAR method [70] to produce an average velocity map by first converting the phase-corrected interferograms into linear displacement velocity maps (mm/year).Then, we calibrate them with an offset value obtained at a stable reference area(s) over one of the velocity maps.ALOS-1 PALSAR Stacking-InSAR is processed separately for the two consecutive ALOS-1 PALSAR frames (see Table I).However, the small overlapping coverage between the two frames does not offer reliable common reference points.Therefore, we choose a single independent area for each ALOS-1 acquisition frame.Then, we average these calibrated LOS displacement velocity maps.The "random" noise due to the atmosphere interferometric phase component fades by averaging.Finally, we convert this LOS product into average vertical velocity, assuming that the horizontal deformation is null.We apply the same stacking method to the ALOS-2 interferograms.To calibrate the ALOS-2 LOS velocity maps, we apply an offset obtained by averaging the two reference areas picked for the ALOS-1 interferogram sets.In Section V-D4, we propose other stacking approaches and their limits.

B. COASTSAT: Detection of Shoreline Retreat
To obtain a time-series of sandy coastline position and understand its evolution, we use the open-source software toolkit COASTSAT developed by Vos et al. [71].COASTSAT retrieves publicly available orthorectified satellite imagery from Google Earth Engine, such as Landsat-5, Landsat-7, Landsat-8, and Sentinel-2, over our region of interest.Table III sums up the multispectral acquisitions used in this shoreline retrieval.We preprocess the dataset by masking the cloudy pixels, applying pansharpening, and downsampling the images.No additional coregistration is applied.We combine a supervised image classification with a subpixel resolution border segmentation to detect the shoreline on each image.We also use COASTSAT to perform the classification using a neural network supervised classifier (i.e., multilayer perceptron from scikit-learn with an error estimates of 1% [71]) to extract the "sand," "water," "white-water," also called surf zone, and "other land features" classes.
From these classified images, segmentation is carried out for the "sand" and "water" classes to locate the coastline.We then apply the modified normalized difference water index (MNDWI), derived from the same multispectral images, and Otsu's thresholding algorithm [72], followed by a contour tracking using a super-resolution enhancement of the marching squares algorithm [71], [73] to improve the precision in locating the coastline.Then, we extract the shoreline time series for the given selected transects (see Fig. 2) after removing the tidal component at the following position: Long, Lat: 50.37177700, 26.17116700.Most of the shores in Bahrain are protected with coastal infrastructure to stabilize coastlines and mitigate beach erosion, e.g., seawalls and jetties.In light of this, from optical satellite imagery observations, we identify two seemingly unaffected beaches, Balkiya and South-West beaches (see Fig. 2), rather than averaging hundreds of transects influenced by such infrastructures.As a result, we choose three representative transects to monitor the shoreline changes from 1985 to 2021.

A. Ground Deformation
Fig. 4 shows the LOS displacement velocities and accelerations derived from SBAS-InSAR with the C-band Sentinel-1A ascending and descending images for the 2016-2019 and 2017-2021 periods (see Table I).The LOS velocity maps show a 7-km-wide by 12-km-long subsidence over the Awali oil field, with a rate of 12-16 mm/year with accelerations up to 6 mm/year 2 .
As mentioned in the methodology, a cubic model would be the best temporal fit for this Sentinel-1A LOS displacement time series.Fig. 5 shows the multitemporal coherence that characterizes the fitting quality, with values ranging from 0 to 1, with 1 being a perfect fit.The multitemporal coherence for the subsiding area ranges between ∼0.3 and 0.6.These values show a medium confidence level for the time-series quality and may suggest a nonlinearity of the ground displacement during this four-year observation (see comments in Section V-D3).
Fig. 6 shows the vertical and eastward displacements extracted from two independent SBAS-InSAR Sentinel-1A analyses with different LOS geometries processed for the same period between March 2017 and December 2019 (see Table II).This vertical velocity is similar in magnitude and spatial pattern to the LOS displacements, as shown in Fig. 4, with a vertical velocity reaching −20 mm/year.Moreover, Fig. 6 shows that the western flank of the subsidence bowl slides eastward with a 3-6 mm/year horizontal displacement and a −3 to −18 mm/year vertical displacement.Furthermore, near point A, the western flank moves eastward, but the eastern flank also moves westward.This horizontal pattern is typical in a subsidence bowl.Fuhrmann and Garthwaite [69] model it with synthetic data.I).The two ALOS-1 PALSAR L-band datasets are processed independently; Frame 510 on the upper half and frame 500 on the lower half.At the center of the subsidence bowl, we measure an apparent LOS deformation of −5 and −18 mm/year for ALOS-1 and ALOS-2 PALSAR, respectively, which correspond to a vertical deformation of −5 mm/year (±1 mm/year) and −21 mm/year (±6 mm/year).Smoke from the chimney generates a false positive at E, corresponding to an apparent uplift of 10 mm/year.Fig. 7 also shows a subsidence pattern over the Awali field for 2007-2010 and 2016-2021, resulting from ALOS-1 PALSAR and ALOS-2 PALSAR ScanSAR L-band Stacking-InSAR analysis, respectively.For the first period, we observe a maximum subsidence of 5 mm/year (±1 mm/year) and 21 mm/year (±6 mm/year) for the second.
From 2016 to 2022, the subsidence pattern and amplitude are similar for both SAR sensors (i.e., Sentinel-1 and ALOS-2 PALSAR), despite their difference in central frequency, period of acquisitions, and deformation retrieval method (i.e., SBAS-InSAR and Stacking-InSAR).These similarities significantly increases the confidence level of the results despite the mediumquality model fitting observed in the SBAS-InSAR.
For the 2007-2010 period, the subsidence bowl, observed with Stacking-InSAR ALOS-1 PALSAR, is 1.5 × 1 km only, and its magnitude is in the same order as the uncertainty error.Therefore, we are uncertain if this latter displacement is real or a residual atmospheric noise.
The subsidence acceleration observed between the 2007-2010 and 2016-2022 periods correlates to the 47% oil production increase between mid-2010 and mid-2011 (see Fig. 1).

B. Shoreline Evolution
Fig. 9 shows the shoreline change at three locations in Bahrain: Malkiya beach (transect #6) and South-West beach #1 and #2 (transects #4 and #2).We compare our shoreline Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.observations extracted with the COASTSAT toolkit for each transect with the ones from [74].
At Transect #2 of South-West beach #2, we observe a shoreline recession of 5 m/year [±9 m root-mean-square error (RMSE)] for 2003-2014 compared with 3.5 m/year (±8.5 m RMSE) reported by Luijendijk et al. [74].However, linear trends could not be established for our other COASTSAT analyses and those conducted by Luijendijk et al. [74] at South-West beach #1 and Malkiya beach due to an inadequate coefficient of determination (R 2 = 0).Nevertheless, we still measure the data span for each transect by calculating the 90% central range R 90c between the 5th and 95th percentiles.Our study and Luijendijk et al. [74] show similar ranges of R 90c for each transect.The shoreline at Malkiya beach remains stable within an R 90c of 9-18 m range.Similarly, the shoreline of South-West beach #1 transect #4 is also stable between 1985 and 2021, with R 90c ranging between 20 and 30 m.At transect #2, the shoreline remains stable from 1985 to 2002 within an R 90c range of 13-22 m.

A. Origins of the Subsidence
Our investigation aims at assessing Bahrain's nationwide subsidence, its drivers, the limitation in our approach, and the suggested path forward to reduce the ambiguities with the above.
Our first hypothesis arises from the active oil production occurring at this specific location, where several factors would create residual subsidence (see details in Section VI-A).The second hypothesis is associated with the depletion of the aquifers.We suggest that the Rus-UER aquifer could be semiconfined locally and not unconfined as mentioned above (see cross section B-B' in Fig. 3).Hence, the decrease in pore pressure in the aquifer would increase the effective stress in the upper confining bed (shales) that could compact.
Our third hypothesis suggests that the anhydrites and chalky limestones in the Rus formation could dissolve in contact with groundwater, creating cavities that would collapse when water levels decline.However, the corresponding InSAR displacement pattern would significantly differ from the one observed in Bahrain, with higher displacement frequencies expected in time and space [75].If the first two hypotheses are invalid, we suggest exploring the third hypothesis to determine the presence of active karstification.

B. Forecasting Subsidence for 2100
The deformation maps derived from InSAR show a 7 × 12 km surface subsidence centered over the Awali oil field.The vertical deformation rate reaches a maximum of 18-20 mm/year at the center of the subsidence bowl, while the coastal areas seem stable.We describe, in Section V-D3, the reason for hypothesizing the nonlinearity of the subsidence when assessing ground displacements.As such, predicting future subsidence is challenging.However, as a first-order approximation, one would expect an additional subsidence of 1.4-1.6 m by 2100 if the deformation is constant and linear up to that year.

C. Impacts of the Subsidence
Coastal communities adapt differently to subsidences, flooding, and SLR.For instance, the adjacent cities of Wilmington and Long Beach (Los Angeles-California area) had engaged in 1958 in a vast program of water injection to stop subsidence Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.associated with oil extraction, which reached a total of ∼9 m [76].On the contrary, Jakarta stakeholders implemented a much different strategy.The latter city is now struggling to mitigate its groundwater-related subsidence, i.e., up to 20-28 cm/year, with concrete dikes and seawalls [12].Despite these efforts, the Indonesian government decided in 2019 to move its seat to East Kalimantan on Borneo Island [12].
Monitoring subsidence is crucial to understanding its amplitude and extent and forecasting flooding vulnerability.On the small island of Bahrain, our subsidence mapping can improve the prediction of future vulnerable lowlands and modeling of local relative SLR and shoreline retreat.
1) On Coastal Flooding: The subsidence observed in Bahrain spans 7 km over an island that is 14 km wide.Low-lying coastal zones are susceptible to seawater flooding, particularly after storm surges (see the 1 m above mean sea-level (±2 m) in Fig. 2).Any subsidence adjacent to these areas can worsen the flooding, resulting in wider and more prolonged inundation periods.
Seawater intrusion arising from aquifer depletion is already an issue in Bahrain.How would the expansion of flooding areas impact groundwater?
2) On Oil Production: Faults, fractures, anticline structure, and permeable layers heterogeneity control the gas cap distribution and water ingress in reservoirs [53].These fractures can develop, and new ones can appear when ground subsidence occurs over oil fields [3].Can this mechanical process intensify the seawater intrusion into the aquifers?

D. InSAR Limits 1) Errors in Assessing Displacements:
We estimate the uncorrelated errors corresponding to the inherent noise of the InSAR data that is comparable at every point.This noise contains the thermal noise and remaining RFI residuals from each SAR acquisition, spatial and temporal decorrelations at surface scatterers, and the soil moisture change [77] between each interferometric pair.
To assess the uncorrelated errors for Sentinel-1A SBAS displacement time series and velocity maps, we first determine the variability in the temporal displacements over a representative "stable" area.We calculate this error estimation by adding the median with half of the 90% central displacement range (R 90c ) between the 5th and 95th percentiles.For the Sentinel-1A time series (see Fig. 10), we choose a "stable" point at N, away from the SBAS-InSAR reference point F (see Fig. 4), where the 2016-2021 median vertical displacement is null and constant (± 1 mm).At N, the median vertical displacement is 2.5 and 3 mm, spanning over 5 mm (R 90c ), for the whole period, respectively, for ascending and descending directions.Using these values, we evaluate a displacement error of ±5 mm (see calculations in Fig. 10).We also measure the error at E, where the smoke arising from the aluminum factory is persistent for each SAR acquisition and affects the InSAR results.At E, the median vertical displacements are 1.5 and 1 mm, R 90c are 10 and 12 mm, and error estimates are ±7 mm for the ascending and descending directions, respectively.Fig. 6 shows the Sentinel-1A vertical and eastward velocities, where the red dashed square represents a "stable" region.Fig. 12 shows the corresponding histograms of the variability in spatial displacements over this "stable" area.The median vertical and eastward velocities are −1 mm/year and 0 mm/year, respectively, and their data range R 90c is 3 mm and 1.5 mm for the vertical and eastward velocities, respectively.Therefore, we determine a vertical velocity error of ±3 mm/year and an eastward error of ±1 mm/year.
To estimate the uncorrelated error of the ALOS-1 PALSAR (2007-2010) and ALOS-2 PALSAR (2016-2021) velocity maps (see Fig. 7), we select pixel values outside the subsidence zone, distributed over the whole map, i.e., 137 and 200 pixels for ALOS-1 PALSAR and ALOS-2 PALSAR, respectively.Fig. 11 shows the corresponding histograms.For these "stable" points, the ALOS-1 PALSAR vertical velocity median is 0 mm/year, and R 90c is 2 mm.Hence, for this 2007-2010 velocity map, we determine a vertical velocity error of ±1 mm/year.The "stable" points for the ALOS-2 PALSAR vertical velocity have a median and R 90c , of 1.5 mm/year and 8 mm, respectively; hence, we determine a vertical velocity error of ±6 mm/year for the 2016-2021 velocity map.
On the contrary, the correlated errors are associated with inaccuracies in the SBAS-InSAR deformation retrieval models.The quality and magnitude of this error are assessed using the model fitting coefficient (or multitemporal coherence) and RMSE (see Fig. 5) between the InSAR observations and the deformation model.The RMSE, in particular, serves as a firstorder estimate of the correlated errors.Ancillary datasets, such as ground deformation control points derived from GPS, are an important support to constrain and validate the InSAR observed displacements.However, the study area does not have a publicly available GPS dataset or continuous geodetic data to support our analysis.The RMSE is the highest at E, where the smoke is persistent (see Fig. 5) and reaches 2-3 mm, whereas it ranges between 0 and 2 mm in other places.
Equation ( 2) expresses the stochastic components of the observed InSAR unwrapped phase [see (3)] [78].Each component is characterized by independent physical processes and, thus, has singular statistic characteristics and errors [78] where ΔФ is the unwrapped phase at pixel x between two SAR acquisitions.The unwrapping processing may cause a residual phase [79], but Agram and Simons [78] suggest that current sensors detect sufficient scatterers at the resolution cell to unwrap the phase accurately.ΔФ DEF is the cumulative surface deformation in the time spanning the SAR acquisitions.ΔФ ATM is the difference in propagation delay through the troposphere and in ionosphere effects between the acquisitions.ΔФ SCATT is the change in surface scattering.ΔФ UNCORR is the combination of all other uncorrelated noise sources of the interferometric phase [78].The latter term includes InSAR processing errors caused by coregistration and interpolation errors [80].In this equation, we omit the phase contribution corresponding to orbit and topographic errors since they are well estimated with precise orbital information [78] Φ where Ф(x) is the estimated absolute phase value (unwrapped), ϕ(x) is the measured modulo-2π phase value (wrapped), k(x) is an integer value, and x is a given pixel [79].
The retrieval of ΔФ is achieved with the InSAR time series, where the following points are noticed.
1) Ф DEF is correlated in space and in time.
2) Ф ATM is correlated in space, showing similar values in adjacent pixels, but is uncorrelated in time because of the "large" temporal resolution of Sentinel-1A, i.e., 12 days.3) Ф SCATT is correlated within the resolution pixel for interferograms having a SAR acquisition in common but is uncorrelated temporally [80] and spatially compared with adjacent pixels [78].4) The RFI residuals are correlated in space, but we are unsure if they are correlated in time.Finally, the ground deformation is likely isotropic, meaning that its northward deformation component does not affect the estimation of the vertical and eastward displacements.Fuhrmann and Garthwaite [69] suggest that northward displacement residuals would only contribute to the inaccuracy of << 1 mm/year on the vertical displacement.
2) Absolute Elevation of the Subsidence: The InSAR-based ground deformation model, referred to as the relative ground deformation, relies on reference points selected outside the subsidence region in an area assumed stable during the study period.To accurately determine the absolute ground deformation allowing comparison with shoreline retreat, the InSAR model must be calibrated with GPS station(s).However, this limitation is not critical for our study as the recent land reclamation appears to prevent Bahrain's shoreline recession.
3) SBAS Model Fitting: Fig. 5 shows the multitemporal coherence, which is a metric of the model fitting with the InSAR observations.As the time series lengthens, the model fitting becomes more challenging.This issue may be attributed to the nonlinear displacement caused by oil pumping, as illustrated in the acceleration displacement maps in Fig. 4. It can also be attributed to smoke residuals resulting from the associated petroleum gas (APG) combustion (see Section V-D5).
To improve the SBAS model fitting and increase the accuracy of time-series displacement, we recommend splitting the dataset into periods of approximately two years and rerun the SBAS analysis.
4) Atmosphere Removal With Stacking-InSAR: Our Stacking-InSAR analysis converts each interferogram into a velocity displacement map and averages them.This method exaggerates the atmospheric effects when using interferograms with small temporal baselines.It also limits the weight of the apparent deformation when using interferograms with large temporal baselines.The performance of this method drops when the deformation is minor, and the radar wavelength is large.
For instance, this method would not be appropriate when combining 1-month and 2-year L-band interferograms to characterize minor subsidence of 5 mm/year.In this case, we could ponder the interferograms as a function of their temporal baseline.
However, this latter method would intensify the atmospheric noise from "noisy" interferograms.
For ALOS-1 PALSAR and ALOS-2 PALSAR Stacking-InSAR, we carefully chose a large temporal baseline for all the interferograms to favor the deformation over the atmospheric noise.
5) Water Vapor Plumes From Industrial Chimneys: We measure that the vertical SBAS-InSAR measure uncertainty (i.e., temporal variability) near an aluminum factory is ±7 mm (at E) with C-band Sentinel-1A, for the 2016-2022 period, in comparison with the overall noise level estimated at ±5 mm (at N).At E, we believe that the smoke continually emanates from the chimney(s) and that the plume follows the preferential wind directions.In fact, the wind in Bahrain typically blows from the North-West throughout the year and from the South and South-East in July and August, with an average speed of 4-5 m/s (σ = 3) [81].If this plume is persistent in time and space, the SBAS algorithm will fit a deformation model to the corresponding InSAR noise.To illustrate this issue, we select two SAR scenes (October 2019 and 2020) from the SBAS time series.As shown in Fig. 13, the plumes in the October 2019 and October 2020 scenes have dimensions of 10 × 3 km and 10 × 5 km, respectively, originating from East-North-East and North-East winds, respectively.The resulting phase delay of 1-2π radians is removed with the SBAS atmosphere removal.The SBAS-InSAR effectively filters out smoke at E, although there is a small residual error with an RMSE of 2-3 mm.However, the Stacking-InSAR LOS displacement velocity maps exhibit persistent noise at E, as shown in Fig. 7.This figure reveals an apparent local uplift of ∼6 mm/year with ALOS-1 PALSAR for 2007-2010 and greater than 10 mm/year with ALOS-2 PALSAR for 2016-2021.
In the Awali oil field, numerous small degassing chimneys are scattered throughout the area.These chimneys emit flare, smoke, and water vapor resulting from the combustion of APG, a waste product of petroleum extraction.These chimneys appear much smaller than the factory's chimney at E, and we are unsure if they are frequently active.Nowadays, to reduce the carbon footprint of oil extraction, APG is often utilized to power local facilities [82].Despite the low likelihood of flares impacting the SBAS-InSAR, we reprocess the Sentinel-1A datasets to address any potential effects of local, high temporal variability below the 12-day Sentinel-1A revisit time.To further address this issue, we implement a spatial filtering of 500 m and a temporal filtering of ten days for the SBAS atmosphere removal to accurately resolve the spatial and temporal distribution of flare/smoke variations.For our initial SBAS atmosphere removal process, we implemented a spatial filtering of 1400 m and a temporal filtering of 600 days to smooth spatially and temporally the InSAR time series.To conclude, both SBAS analysis shows the same subsidence pattern in extent and amplitude, which confirm that our first SBAS was relevant.
We recommend identifying plumes and masking the affected pixels on interferograms by comparing the SAR acquisitions with multispectral imagery that was acquired simultaneously.Next, one can conduct the intermittent SBAS-InSAR analysis, which integrate these intermittent coherent targets, affected by smoke, in the deformation retrieval.

E. Shoreline Change
In the Gulf, shoreline fluctuates in response to sea-level changes and coastal sediment transport dynamics, which occur naturally but can be accelerated by anthropogenic activities, such as dredging, coastal defense, and land reclamation, resulting in either beach accretion or erosion.
Aladwani [83] measured shoreline changes in the southern shore of Kuwait based on shoreline detection analysis similar to the one performed in this study with the COASTSAT toolkit.Their study reveals a succession of accretion and erosion from 1986 to 2021, i.e., accretion of 1.13 m/year for 1986-1995, erosion of 0.20 m/year for 1995-2005 and 1.31 m/year for 2005-2015, and accretion of 1.49 m/year for 2015-2021.Similarly, Zeinali et al. [84] monitored the shoreline changes for the Boushehr Province in the southern part of Iran in the Persian Gulf.Their study covers 64 km of coastline and includes 642 transects.It reveals that erosion affects 60% of the coast, with average rates of 0.3 m/year in Shirino, 0.5 m/year in Nayband Bay, and 0.9 m/year in Haleh.Moreover, they observed accretion with average rates of 17.8 m/year in Nakhl Taghi and 0.4 m/year in Asaluyeh.We note that neither study refers to SLR when describing shoreline changes but rather focuses on erosion/accretion.
Characterizing shoreline change solely in the context of SLR is inadequate and requires a comprehensive understanding of the underlying sedimentary processes.Therefore, to ensure the accurate representation of the shoreline in natural conditions and deconvolve the contribution of the SLR on shoreline dynamics, we select three COASTSAT transects that remain unaffected by coastal defense or land reclamation activities.We then compare these transects over the 1985-2022 period with those monitored by Luijendijk et al. [74] (see Fig. 9).Averaging shoreline changes from hundreds of transects, such as Aladwani [83], in the southern shore of Kuwait and Zeinali et al. [84] in the Boushehr Province of Iran, would not be a representative approach in the case of Bahrain.
Stable shorelines are observed at the three selected locations, except for Transect #2 of South-West beach #2, which shows a shoreline retreat of 5 m/year (±9 m RMSE) for 2003-2014.This high rate of shoreline retreat is consistent with findings reported by Luijendijk et al. [74], who observed a rate of 3.5 m/year (±8.5 m RMSE) for the same location and time period.This rate is significantly higher than the average shoreline retreat of ∼0.80 m/year (see details in Appendix B) derived from the local SLR obtained from tide gauge measures and corrected for VLM.
Further investigation is required to determine whether this 11-year-long shoreline retreat, occurring between 2003 and 2014, results from sediment transport dynamics or relative SLR, which may be associated with the nationwide subsidence centered over the Awali oil field that started between 2010 and 2016.

VI. RECOMMENDATIONS AND FUTURE PROSPECTS
We dedicate this section to suggest future work in Bahrain regarding the oil field and the possible ground deformation in the northern part.

A. Geologic-Hydrogeologic Deformation Model
Kuzmin [85] states that apparent subsidence over an oil field could be a residual deformation involving sequences of uplifts and subsidences.This statement agrees with our hypothesis on a nonlinear deformation over the Awali field (see Section V-D-3).
To model this residual deformation (i.e., the sum of the deformations) over the Awali field, one would need to take into consideration the following: 1) volume of crude oil extracted in each reservoir and change (decrease) in the reservoir pressure [85]; 2) dynamics (decrease) of poroelastic parameters and compressibility of the rocks [85]; 3) weight of the overlying rocks [85]; 4) tectonic (anticlinal uplift) [85]; 5) volume of gas or water injected to maintain the pressure and large gas cap to increase the oil production [52]; 6) pore pressure change of the aquifers; 7) variation of the hydraulic head [86]; 8) volume of water encroaching.However, before performing such complex modeling, one could start by comparing our InSAR deformation observations with a first-order deformation approximation using the Mogi model [69] of the main oil reservoir, i.e., the Mauddud reservoir.
The accurate structural and hydrogeologic deformation model would provide unique insights to address the following questions: First, does the deep compaction from the oil reservoir affect the shallow groundwater storage?Second, does the subsidence increase the seawater intrusion into aquifers?

B. InSAR Over the North of Bahrain
As explained in Section V, the prediction of future lowlands and shoreline retreats can be improved by monitoring the subsidence in Bahrain.The northern part is the most inhabited part of the country, and this region should also be investigated to complete our work on local relative SLR implications.In 2005, Al Zubari [54] reported that the confined Dammam aquifer water levels were declining in the North of Bahrain while the salinity was increasing.Therefore, we suggest sensing possible subsidence arising from the aquifer depletion and consolidation of the aquitards.
However, we stress some limits in performing InSAR in this area.First, during the "InSAR reflattening" processing step, we focus on the central part of Bahrain; thus, ground deformation could have been underestimated in northern Bahrain, i.e., from the ALOS-1 and ALOS-2 velocity maps, we do not observe ground deformation in that part of Bahrain.We suggest to subset the single look complexes (SLCs) over the North before producing the interferograms.Second, we warn that the RFI was strong over that area in 2016-2022; thus, it will be difficult and time-consuming to remove it to produce InSAR time series (see Fig. 14 in the Appendix E).However, RFI may be absent or rare in the 1990s-2000s.We suggest using older SAR sensors, such as ERS-1, ERS-2, ENVISAT, and Radarsat-1.Third, spatial filtering should be executed with vigilance.In fact, InSAR measures the average ground displacement with the main lobe of the radar antenna, sensing all the scatterers located in the resolution cell.In cities, several strong scatterers can occur within the resolution cell.Hence, the scattering arising from sidelobes interferes with the adjacent scattering of the main lobe.Despite their long computation time, the permanent scatterers InSAR techniques [87] solve this issue by performing phase interferometry for each scatterer, then selecting only the ones sensed by the main lobe.On the other hand, DinSAR techniques can be challenging and deserve attention.For instance, for this particular issue, the original Goldstein adaptive filter produces cross-like patterns that limit the use of InSAR when not tuned with care [88].

VII. CONCLUSION
Between mid-2010 and mid-2011, crude oil production in Bahrain increased by 47%.This development prompted our investigation into ground deformation over the Awali oil field, as the fluid withdrawal may result in reservoir compaction and land subsidence.In addition, our study aims to explore the potential implications on relative SLR and shoreline retreat.
In this study, we use the Stacking-InSAR method with L-band ALOS-1 PALSAR and ALOS-2 PALSAR ScanSAR SAR satellites to measure the expanse and magnitude of surface deformation in the Awali oil field, Bahrain, from 2007 to 2010.Moreover, we apply the SBAS-InSAR method with the C-band Sentinel-1A SAR satellite in ascending and descending directions to measure deformation from 2016 to 2021.
We observe from both methods a 7 × 12 km subsidence over the Awali oil field with an uncalibrated vertical velocity of up to 20 mm/year between 2016 and 2021.However, we do not observe significant subsidence for 2007-2010.We stress here the possible consequences of 7-km-wide subsidence over the 14-km-wide island of Bahrain.The nationwide subsidence can cause significant seawater intrusion resulting in groundwater contamination and can trigger relative SLR and shoreline retreat.
We recommend performing a similar InSAR analysis in other low-lying Gulf islands where oil exploitation occurs, and relative SLR is a concern when forecasting coastal infrastructure planning.
Despite the current land reclamation policy impeding most of Bahrain's shorelines' recession, we observe from subpixel shoreline detection analysis, a shoreline retreat of up to 5 m/year between 2003 and 2014 in the South-West part of Bahrain.Further investigations are required to determine whether this 2003-2014 retreat results from natural sediment transport dynamics or relative SLR originating from the nationwide subsidence that started between 2010 and 2016.

A. Absolute Sea-Level Rate Estimation
The Gulf region is affected by an atmospheric pressure change of 10 mbar/year, called the inverted barometer (IB) effect, which corresponds to an annual sea-level amplitude variation of 10 cm/year [38].This latter effect causes an ocean bottom deformation variation called ATMospheric Loading (ATML) [38].Eventually, IB and ATML effects cause mass redistribution in the land-atmosphere-ocean system leading to a third variation of the sea-level [38].This latter effect is called self-attraction and loading (SAL) [38].The IB effect, corrected from ATML and SAL, should be subtracted from tide gauge observations, except in Bahrain, where it seems uncorrelated to the sea-level [38].After removing the 18.61-year nodal cycle, Alothman et al. [38] estimate a relative SLR rate of 2.97(±0.51)mm/year at Manama, North-East of Bahrain island for the period 1979-2007 [38].From this rate, we remove the postglacial rebound that is around −0.31 mm/year in this region [89] and obtain the absolute SLR rate of 3.28(±0.51)mm/year.Nevertheless, this rate should be corrected by the VLM.The latter shows a velocity ranging from −0.33(± 0.20) mm/year [39] to −0.51(± 0.24) mm/year between mid-1996 and mid-2008, over the GPS station "BAHR" that is 3 km away from the tide gauge [38].Considering these VLM rates, we estimate a corrected-absolute SLR rate ranging from 2.77(±0.56)mm/year to 2.95(±0.55)mm/year.

C. Uncorrelated Error
We define the uncorrelated error as the sum of all the noises from the SAR data.

D. Effect of Smoke on Phase
In Fig. 13, the plumes with a size of ∼10 × 4 km introduce a 1-2π radians phase delay that can be seen in the interferogram.However, the SBAS atmosphere removal can remove this noise effectively.

E. RFI Removal
Military, airplane, vessels, or weather radar emissions interfere with the SAR radar satellites emitting at a central frequency of 5.405 GHz and 1.27 GHz, respectively, for Sentinel-1A and ALOS-1.It results in RFIs affecting the Level-1 SLC radar observations and the InSAR products.To eliminate this noise, we mask the RFI on the defocused SLCs with a notch filter in this Azimuth-time range-frequency domain, as described in [66]; we then refocus the SLCs.Fig. 14 shows an example of SAR SLC amplitude images affected by RFI.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 10.Vertical displacement time series at N, E, and F. The vertical displacement is estimated from the LOS and the incidence angle.We choose N as quality control over a region where the velocity is null.E is located near the persistent smoke arising from the aluminum factory.F is the relative displacement reference used during the SBAS-InSAR; no displacement should be observed at this location.We calculate a conservative uncertainty error of ±5 mm from the time series at N for both the ascending and descending acquisitions.At E, the error estimate is ±7 mm, and at F, the displacement is null for the whole time series (±1 mm).

F. Ionospheric Noise Removal
Ionospheric noise on SAR images, associated with the Appleton anomaly [90], is expected on scenes acquired during the day at low latitudes, i.e., between −15 and +15° [65].The atmospheric noise removal in the SBAS is not always appropriate to remove such spatial high-frequency noise.Fig. 15 shows Fig. 11.Histograms of the vertical and eastward displacements variabilities over a stable area (see Fig. 6).From 5th and 95th percentiles, we measure a vertical error ranging from mm/year 2.2 to +1.0 mm/year and an eastward error ranging from −0.9 to +0.6 mm/year.We determine a conservative vertical velocity error of ±3 mm/year and an eastward error of ±1 mm/year for the Sentinel-1A LOS combination.Fig. 12. Histogram of the vertical displacement over stable areas for ALOS-1 PALSAR and ALOS-2 PALSAR Stacking-InSAR velocity maps (estimated from the LOS displacement and the incidence angle).For ALOS-1, the median is 0 mm/year, and the 5th and 95th percentile are −1 mm/year and 1 mm/year, respectively.For ALOS-2, the median is 1.5 mm/year, and the 5th and 95th percentile are −2.5 mm/year and 5.5 mm/year, respectively.The conservative vertical velocity errors are ±1 mm/year and ±6 mm/year for ALOS-1 and ALOS-2 PALSAR, respectively.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.sharp phase changes related to the ionosphere phase difference, which can create phase islands and unwrapping errors during the unwrapping sequence of the InSAR processing [91].We propose the following approach to solve this issue.First, mask the areas of rapid phase change on the affected interferograms and then apply intermittent SBAS-InSAR [67], which computes temporal interpolation between coherent pixels and estimates the ground displacement over these masked areas.
We did not have to follow this approach in this work since only a few interferograms were affected, for which the SBAS atmospheric noise removal was effective.

Fig. 1 .
Fig. 1. Volume of crude oil produced in Bahrain, mainly from the Awali oil field, between January 2009 and December 2019 (open database) [45].The production had increased by 47% between mid-2010 and mid-2011.

Fig. 2 .
Fig. 2. Geologic map of Bahrain overlaying the island DEM (top frame) and topographic profile of the Awali oil field in the central depression of Bahrain (bottom frame).

Fig. 5 .
Fig. 5. (a) and (c) Multitemporal coherence and the (b) and (d) RMSE maps allow us to control the model fitting quality of the four years of Sentinel-1A displacement velocities.The quality of the cubic fitting is medium over the subsidence bowl, with a multitemporal coherence ranging from 0.35 to 0.5.The RMSE is adequate, with an error under 2 mm.

Fig. 7 .
Fig. 7. Stacking-InSAR LOS displacement velocities for 2007-2010 and 2016-2021 from the sensors ALOS-1 and ALOS-2 PALSAR L-band, respectively (see TableI).The two ALOS-1 PALSAR L-band datasets are processed independently; Frame 510 on the upper half and frame 500 on the lower half.At the center of the subsidence bowl, we measure an apparent LOS deformation of −5 and −18 mm/year for ALOS-1 and ALOS-2 PALSAR, respectively, which correspond to a vertical deformation of −5 mm/year (±1 mm/year) and −21 mm/year (±6 mm/year).Smoke from the chimney generates a false positive at E, corresponding to an apparent uplift of 10 mm/year.

Fig. 9 .
Fig. 9. Shoreline change occurring at the three locations Malkiya beach (frames a and b), South-West beach #1 (frames c and d), and South-West beach #2 (frames e and f) (see locations in Fig. 2).For each location, we compare our shoreline change observations in frames (a), (c), and (e), with those of [74] in frames (b), (d), and (f).At Malkiya beach and South-West beach #1, we cannot express a shoreline change rate with linear regression (R 2 = 0); instead, we give the 90% central range R 90c representing the amplitude between the 5th and 95th percentile.Therefore, at these two locations, we observe for 1985-2021 that the shoreline is stable with an uncertainty of R 90c = ∼20 m.At South-West beach #2, the shoreline is also stable for 1985-2002 (R 90c = ∼22 m) and 2015-2021 (R 90c = ∼34 m); however, we observe a shoreline retreat for 2003-2014 with a rate up to 5 m/year (RMSE ±9 m).

Fig. 13 .
Fig. 13.Phase delay from smoke (left) and removal during the SBAS (middle and right).

Fig. 14 .
Fig. 14.RFI removal.The 71 Sentinel-1A scenes are affected in this study and need RFI removal.(a) Scene not affected.(b)-(d) Scenes affected by RFI.

Fig. 15 .
Fig. 15.Example of ionospheric noise removal.On the left frame, the 12day interferogram shows a solid ionospheric contribution from the 2017-11-04 acquisition.This noise is well removed using the SBAS-InSAR methodology.