MONITORING LAND SUBSIDENCE IN EGYPT’S NORTHERN WEST COAST USING INTERFEROMETRIC SYNTHETIC APERTURE RADAR.

: The northern west coast of Egypt has experienced significant development and urbanization in response to its growing population. This has resulted in the construction of many new cities along the Mediterranean Sea shore, making the region a major attraction for both economic and recreational purposes. However, this rapid growth has also created a need for continuous monitoring to ensure environmental and economic sustainability and to identify potential hazards. Remote sensing techniques, such as satellite imagery and interferometric synthetic aperture radar (InSAR), have proven to be effective tools for monitoring remote regions. However, InSAR has some limitations, such as susceptibility to tropospheric and ionospheric delays as well as baseline errors. To address these limitations, this study utilized a differential technique that uses the phase difference between two nearby pixels to reduce the spatially correlated atmospheric and baseline errors. In this study, Sentinel-1 SAR images were used to create a single master InSAR stack covering the study region from 2018 to 2023. The estimated land surface velocities showed a subsidence motion of 10 to 15 mm/year in the region of Marsa Matruh city, which is consistent with existing urbanization and agricultural activities. The findings of this study highlight the importance of continuous monitoring and the potential of remote sensing techniques for identifying potential hazards and ensuring sustainable development.


INTRODUCTION
In recent years, Egypt has undergone significant development and urbanization to address the needs of its growing population.Many new cities have been constructed throughout the country, with one of the largest regions being the northern west coast of Egypt.This region has witnessed the construction of numerous cities along the Mediterranean Sea shore, extending from Alexandria city to Marsa Matruh city, including New Alamin city (Attia et al., 2019).Additionally, the ElDabaa nuclear plant is being established in this region (Megahed, 2009).Large areas of land have also been prepared for agriculture, and new roads and railroads have been built to meet the transportation demands of the population.As a result of these developments, this once-desert region of Egypt has become a major attraction for many citizens due to its economic and recreational offerings.However, the rapid urbanization of the region also emphasizes the need for continuous monitoring to ensure environmental and economic sustainability as well as identify potential hazards that may arise due to this development, (Iskander, 2021).Therefore, it is crucial to maintain the environmental and economic sustainability of the region by undertaking continuous monitoring and identifying any potential hazards that may arise from this rapid urbanization.
Satellite remote sensing technology offers a range of applications and solutions for the ongoing monitoring of remote areas.Optical remote sensing, for example, can be used to monitor urbanization growth (ElGharbawi et al., 2022), vegetation (ElGharbawi et al., 2023), and the impact of urbanization on land surface temperature (Moazzam et al., 2022).Another remote sensing method, Interferometric Synthetic Aperture Radar (InSAR), has proven effective in analyzing surface deformation caused by various phenomena, including earthquakes (Massonnet et al., 1994) (ElGharbawi andTamura, 2014) (ElGharbawi and Tamura, 2015), glacier movements (Goldstein et al., 1993), land subsidence (Buckley et al., 2003), and urban hazard as-sessment (ElGharbawi, 2019).However, InSAR has its limitations.One limitation is the susceptibility of InSAR observations to tropospheric and ionospheric delays, which can significantly impact accuracy.The quality of satellite orbital data also has a significant effect on InSAR observations.In addition, deformation values obtained from InSAR can be affected by phase unwrapping errors and present relative deformation values.Despite these limitations, InSAR remains a valuable tool for monitoring remote areas and is widely utilized in various applications.
Several approaches have been proposed to reduce the tropospheric delay artifact observed in InSAR.For example, (Li et al., 2006) and (Song et al., 2008) combined a topographydependent turbulence model with the ZTD derived from GPS, while (Onn and Zebker, 2006) used an exponential altitudedependent model derived from GPS ZTD.In the Tokyo area, Japan, (ElGharbawi and Tamura, 2014) and (Janssen et al., 2004) used the continuous GPS network's ZTD and the kriging algorithm to correct for the tropospheric delay.
For baseline errors, (Zhang et al., 2008) applied a nonlinear algorithm to unwrapped phase maps and a bilinear approach to interferograms.In the southern San Andreas fault, (Fialko, 2006) removed a linear ramp from a stack of interferograms using GPS velocities.These techniques have proven effective in reducing the impact of tropospheric and baseline errors on In-SAR readings, improving the accuracy of surface deformation measurements.
Since there is no permanent GPS station present in our study region, it is not possible to carry out a full GPS-based atmospheric and baseline correction.To address this, we utilized a differential technique that took advantage of the phase difference between two nearby pixels to reduce spatially correlated atmospheric and baseline errors.This technique was originally proposed for InSAR permanent scatterers by (Ferretti et al., 2000), and it has been successfully applied for estimating land surface deformation in vegetated and rough terrain regions using InSAR distributed scatterers (ElGharbawi and Tamura, 2022).By utilizing this differential technique, we were able to effectively reduce the impact of atmospheric and baseline errors on our InSAR observations, improving the accuracy of our measurements.
This study utilized InSAR technology and Sentinel-1 SAR images to monitor the mean velocity of land deformation along Egypt's western north shore over a period of five years.We created a single master InSAR stack covering the study region from 2018 to 2023 using Sentinel-1 SAR images.In the absence of GPS stations in our study area, we used the Generic Atmospheric Correction Online Service for InSAR (GACOS) (Yu et al., 2018) as a first correction for the tropospheric artifact in the generated interferograms.Since Sentinel-1 employs a Cband signal that is minimally impacted by ionospheric activity, we did not address the ionospheric artifact.
The estimated land surface velocities reveal a subsidence motion in the region of Marsa Matruh city that is consistent with the existing urbanization and agricultural activities in the area.This finding highlights the impact of human activities on land deformation, emphasizing the importance of continuous monitoring to ensure sustainable development in the region.

Study Area
This paper focuses on the Northwestern Mediterranean Coast of Egypt (Figure 1), which has undergone significant urbanization in recent years.This coastal region, which stretches from 30.5 o to 31.5 o North and from 26.5 o to 29 o East along the Mediterranean Sea (Figure 1.a), features a unique blend of urban and natural elements.The coastline is characterized by lagoons, rocky cliffs, and sandy beaches that serve as habitats for various marine species.Additionally, the area contains several wetlands that are home to migratory birds and local fishing populations.
On the other hand, the urbanization of the northern west coast of Egypt is characterized by towering apartment buildings, sprawling resorts, and vast industrial facilities.This contrast between urbanization and natural features highlights the need for continuous monitoring to ensure sustainable development and to identify potential hazards in the region.

Dataset
To conduct our analysis, we utilized data from Copernicus Sentinel-1, processed by the European Space Agency (ESA).Specifically, we used the vertical-vertical polarization (VV) mode of Sentinel 1-A's C-band Interferometric Wide swath (IW) SAR images to gather eleven pictures covering the study area (Figure 1.b).We employed single-look complex (SLC) images with multi-looking of 4 by 20 pixels in both the azimuth and range directions.
The final results were geocoded in Universal Transverse Mercator (UTM) zone 35 N using the Copernicus Digital Elevation Model (DEM) 30m resolution (GLO-30) [https : //spacedata.copernicus.eu/].We produced ten single master interferograms and related coherence maps using the acquired images (Table 1).This approach allowed us to perform surface monitoring over a period of five years while maintaining the phase coherence necessary to obtain more accurate estimations. No

METHODOLOGY
This section provides a detailed description of the proposed methodology, which is illustrated in Figure 2. The following steps were taken: (1) We created a single master interferometric stack using the existing SAR acquisitions.(2) To account for atmospheric delay, we utilized GACOS data for preliminary adjustment.
(3) To correct the baseline error for each interferogram, we performed spatial filtering using a 2D surface of the second degree.(4) Using a Delaunay triangulation network, we selected a set of sparse, highly coherent pixels, connected them, and estimated the deformation velocity and DEM error for each arc.We then inverted the solution to estimate the deformation velocity and DEM error for each coherent pixel.(5) We further adjusted the InSAR stack by employing the residual phase of the chosen coherent pixels.( 6) Finally, we selected all feasible pixels in the InSAR stack while estimating their deformation velocities.This methodology enabled us to effectively reduce atmospheric and baseline errors, leading to more accurate estimations of deformation velocities and DEM errors.

InSAR stack generation
Equation 1 presents the interferometric phase resulting from the complex conjugation of two single look complex (SLC) images.
The phase is impacted by topography ϕT opo, which is removed using DEM.In this study, we utilized the Copernicus GLO-30 DEM to eliminate the topographical impact observed in the interferograms.However, the interferometric phase is also affected by deformation ϕ Def orm , tropospheric delay ϕT ropo, ionospheric delay ϕIono, baseline errors ϕ Baseline , and noise ϕNoise.These effects must be accounted for to obtain accurate measurements of deformation velocities and DEM errors.
Tropospheric delay can vary significantly and is affected by local altitude, temperature, water vapor, and wind.This variability can impact the entire deformation map and produce local disruptions.However, the impact of ionospheric delay can be discounted when using high-frequency radio waves like the Cband used in Sentinel-1, as it is primarily determined by solar activity and the inverse squared frequency of the radio wave.
In Step A of our methodology, we generated complex interferograms (Table 1) (Figure 2 Step A), which were flattened using DEM, filtered, and phase unwrapped using Minimum Cost Flow (MCF).For each interferogram, the unwrapped phase was transformed into displacements using Equation 2.
Where λ is the wavelength and ϕ is the unwrapped phase.All InSAR products were processed by ASF DAAC HyP3 2023 using the hyp3 gamma plugin version 5.2.0 running GAMMA release 20210701 [https : //search.asf.alaska.edu].

Coherence Estimation
Phase coherence is a measure of how comparable two acquisitions' interferometric SAR returns are.The amount of measured coherence reflects how much the target has changed between the acquisition dates of the SAR images.Phase coherence is a measure of the similarity of the interferometric SAR returns from two acquisitions.The coherence measurement reflects how much the target has changed between the capture dates of the SAR images.The complex phase coherence γ is calculated using Equation 3 (Zebker et al., 1992) (Zebker and Chen, 2005).The resulting coherence value ranges from 0 to 1, with 1 indicating perfect coherence and 0 indicating no coherence.Coherence values closer to 1 represent stable targets and are preferred for accurate deformation measurements.
where E ⟨.⟩ is the mathematical expectation, * is the complex conjugate operator, s1 and s2 are complex signals of coregistered SAR images, and |.| is the absolute operator.

Tropospheric Correction
The main objective of tropospheric correction is to compute the differential tropospheric delay map for the unwrapped phase map of interferometric synthetic aperture radar (InSAR).To estimate the tropospheric delay at the observation times of the master and slave images for each interferometric pair, we utilized the Generic Atmospheric Correction Online Service for In-SAR (GACOS).By subtracting the two maps (Equation 4), we can estimate the tropospheric delay and use it to correct the generated interferograms, as shown in (Figure 2 Step B).
To distinguish between stratified and turbulent signals from tropospheric total delays and to produce high spatial resolution zenith total delay maps for use in InSAR correction and other applications, GACOS utilizes the Iterative Tropospheric Decomposition (ITD) model.GACOS provides several crucial attributes, including global accessibility, near-real-time functionality, and simplicity of implementation.Because of the large area of our study region, we use the GACOS data to reduce the effect of tropospheric delay in the InSAR stack, which allows us to use longer arcs for differential analysis because it can be assumed that the atmospheric delay is constant within a spatial range of 1 km (ElGharbawi and Tamura, 2014), and due to the sparsity of the coherent pixels in our study region, we need to use arcs with lengths longer than 1 km in order to achieve full connectivity of the study region.

Baseline Correction
We selected highly coherent pixels from the phase unwrapped InSAR stack and used their LOS unwrapped observations to create a 2D surface of the second degree to represent the baseline error of the area.As we only created a single master InSAR stack, the baseline error is uncorrelated between interferograms, so the model fitting was done for each interferogram.After the initial model fitting, we calculated the fitting error for the pixels, deleted any pixels with errors greater than (1.96 σ), which is equal to the 95% confidence level, and then re-fit the model with the remaining pixels.
To apply the baseline correction, we subtracted this 2D surface from the InSAR unwrapped phase stack, as shown in (Figure 2 Step C).This correction helped to reduce the impact of baseline errors on the InSAR measurements and improve the accuracy of the deformation velocity and DEM error estimations.

Residual Atmospheric Correction
Up to this step, we have corrected a baseline error and made an initial estimate of the atmospheric error, leaving the unwrapped phase stack with residual atmospheric errors.These errors have a spatial correlation but no temporal correlation.To address this, we suggest selecting a sparse, highly coherent pixel that covers the entire study region and connecting it to others using a Delaunay triangulation with an arc length of no more than 5 km.Ideally, the arc length should not exceed 1 km, but we had to choose this length arbitrarily due to the sparsity of pixels in some areas and the need for a fully connected network.Additionally, because we conducted an initial atmospheric correction (Figure 2 Step B), we were able to relax the arc length constraint.
To estimate the mean velocity and DEM error for each coherent pixel, (Figure 2 Step D), let the number of differential interferograms is M, number of arcs is S , number of coherent pixels is K, number of unknowns for each pixel equals 2, identified by pixels mean velocity (V) and DEM error (δh).Also, let the phase difference between two coherent pixels, ∆ϕ i a,b , in the i th unwrapped phase map, (Equation 5).For neighboring pixels, the spatially correlated errors, ∆ϕ For each arc, we use L2-norm to solve for the unknown deformation velocities ∆V and DEM error ∆h by minimizing the error function E, (Equation 6), where j = 1:S. (6) Where ∆Ti is the time difference in years between the master and slave images in interferogram i, B⊥, i is the normal baseline for interferogram i, rj and θj are the range distance and look angle for arc j, respectively.
Finally, the deformation velocity and DEM error for each coherent pixel can be estimated using the least squares inversion, (Equation 7), this inversion is applied for pixel velocities and for DEM error, i.e., two times.Here, [A] is the design matrix relating the arcs (S), with the coherent pixels (K), [L] contains the estimated velocities (DEM error) for all the arcs, and [V] is the unknown velocity (DEM error) for the coherent pixels.It should be noted that (Equation 7) is a large sparse linear system that can be solved by an iterative method based on a bidiagonalization procedure that appears in Matlab as a built-in LSQR function (Paige and Saunders, 1982).

Final Velocity Map Estimation
In Figure 2 Step D, we estimated the mean velocity and DEM error for each coherent pixel of the Delauny Triangulation network.Next, we estimated the residual interferometric phase for each coherent pixel and assumed that it represented the residual atmospheric and baseline errors that could not be removed in the previous steps.We then interpolated these residuals to cover the entire study region for each unwrapped phase and applied the correction to each unwrapped phase, allowing us to estimate the final corrected unwrapped phase maps.
Then we used the final corrected unwrapped phase maps to estimate the mean velocity and DEM error for each pixel directly by minimizing the error function E1 for each pixel (Figure 2 Step E) (Equation 8).

ANALYSIS AND RESULTS
We created an InSAR stack using the Sentinel-1 images shown in (Table 1) using Alaska Satellite Facility Vertix service.The interferograms were flattened using Copernicus GLO-30 DEM, unwrapped using MCF, and geocoded with a spatial resolution of 80 m by 80 m. Figure 3 shows the InSAR stack normal and temporal baselines.We estimated the GACOS tropospheric delay at the time of SAR image acquisition and estimated the tropospheric delay difference using (Equation 4).An example of the atmospheric correction for the first interferogram is presented in Figure 4.a.
The estimated tropospheric delay was used to correct the unwrapped interferometric phase after converting the phase to meters using (Equation 2).
Next, we estimated the minimum coherence map (Figure 4.b) for the entire InSAR stack to isolate the permanent coherent pixels.We then selected a sparse number of pixels covering the entire study region and fitted a 2D surface of the 2nd degree to represent the baseline error for each unwrapped phase.An example of the baseline correction for the first interferogram is presented in Figure 4.c.After correcting the unwrapped phase maps for the baseline error, we selected sparse coherent pixels with coherence values higher than 0.7 for the entire InSAR stack.Then, we connected these pixels using a Delaunay triangulation network (Figure 4.d) and estimated the differential phases between arcs to solve for the coherent pixel velocities and DEM errors using (Equation 5), (Equation 6), and (Equation 7).Then, we estimated the phase residuals for the coherent pixels and interpolated them to cover the entire study region, which we assumed to be the residual atmospheric and baseline errors.An example of the residual atmospheric and baseline errors for the first interferogram is presented in Figure 4.e.
After applying all the correction steps, we used the corrected unwrapped phase stack to estimate the mean velocity of the entire study region using (Equation 8).The final LOS velocity map is presented in Figure 5.a , which shows small variation values consistent with sand transportation in the desert region.
To identify vulnerable regions, we set a velocity threshold of 10 mm/year and defined regions with velocity values exceeding this threshold (Figure 5.b).We found that most of Marsa Matruh city shows a mean velocity of 10 : 15 mm/year for the past 4.7 years, which accumulates to an average of 5.8 cm land subsidence during the study timeframe (Figure 5.c).This subsidence rate can be attributed to the high urbanization rate in the Marsa Matruh region, in addition to the excessive use of groundwater for irrigation purposes, leading to groundwater depletion and land subsidence.

CONCLUSION
The rapid development and urbanization of Egypt's northern west coast, from Alexandria to Marsa Matruh, have brought both economic and recreational opportunities to the region.However, it is crucial to continuously monitor the area to ensure environmental and economic sustainability and to estimate potential hazards due to rapid urbanization.In this study, we used the Interferometric Synthetic Aperture Radar (InSAR) technique to monitor land deformation mean velocity of the western north coast of Egypt for the last five years.The lack of permanent GPS stations in the study region made it challenging to conduct a thorough GPS-based atmospheric and baseline correction.Therefore, a differential technique was proposed and utilized to reduce spatially correlated atmospheric and baseline errors.The results revealed a subsidence motion of 10-15 mm/year from 2018 to 2023 in the region of Marsa Matruh city, which is consistent with the existing urbanization and agricultural activities of the region.This subsidence rate can be attributed to the high urbanization rate in Marsa Matruh, in addition to the excessive use of groundwater for irrigation purposes, leading to groundwater depletion and the presented land subsidence rate.This study highlights the importance of continuous monitoring of urbanization and development areas to ensure their sustainability and mitigate potential hazards.The InSAR technique can be an effective tool for such monitoring, especially in areas where permanent GPS stations are not available.

Figure 1 .
Figure 1.(a) General location of Study area, (b) Study areashowing the SAR images boundary and direction.

Figure 4 .
Figure 4. (a) GACOS atmospheric correction for first interferogram, (b) Estimated minimum coherence along InSAR stack, (c) Estimated baseline error for first interferogram , (d) delauny triangulation network and the selected coherent pixels, (e) final estimated atmosheric correction.

Figure 5 .
Figure 5. (a) Final estimated LOS velocities, (b) regions with LOS velocities exceeding 1.5 cm/year, (c) Regions with subcidence velocityu more than 1.5 cm/year inMarsa Matruh, (d) satellite image of part of the affected region.

Table 1 .
Details of SAR images.