Multi-Sensor Satellite Images for Detecting the Effects of Land-Use Changes on the Archaeological Area of Giza Necropolis, Egypt

: The World Heritage Committee has been meeting to discuss the arrangements of existing World Heritage Sites, and, on 22–26 October, the area from Giza to the Dahshur


Introduction
Cultural heritage sites are considered a wealth of human civilization beyond their value and are rare and irreplaceable places that show the lifestyles of ancient people [1,2].UNESCO defined cultural heritage sites as "the result of life and inspiration, which should be protected and preserved" [3].Cultural heritage regions, encompassing the archaeological and historical built environment and movable heritage, are at risk and still suffer from a lack of specific and systematic risk [4][5][6][7][8][9].In 2005, a study prepared by ICOMOS considered the threats to 641 cultural and mixed world heritage properties due to development projects in addition to the lack of quantitative data on trends between 1994 and 2004 [10].Indeed, the risks that threaten heritage sites include natural (e.g., fires, earthquakes, flooding, tsunami, land and mudslides, winds, tropical storms, land subsidence, etc.) and human factors (e.g., urban development, land drainage, agricultural sprawling, etc.) [11,12].Egypt has multiple cultural sites, objects, and artifacts, as well as historic rock art sites [13].For instance, the pyramid fields from Giza to the Dahshur were accepted as one of the World Heritage Sites in October 1979 [14].Unfortunately, thousands of known and unknown archaeological sites in Egypt are at risk of destruction from urban sprawl and expanding development [15].Specifically, uncontrolled and unplanned urban growth threatens cultural and ecological heritage and are among the major causes of archaeological site loss and damage [16].Urban sprawl is generally regarded as one of the major problems accompanying contemporary urbanization in less-developed countries [17].Around the world, the impacts of urban sprawl include the loss of community and environmental assets, such as cultural heritage sites [18,19], as a result of difficult ground and geo-environmental conditions [20][21][22].
Recently, UNESCO became able to monitor the cultural heritage status via remote sensing data (e.g., satellite imagery) [23,24].Indeed, satellite remote sensing is playing an increasingly important role in the risk detection, uncovering, and documentation of natural and human heritage sites [25].Recently, in addition to archaeological site documentation, archaeological remote sensing and associated technologies (e.g., GIS tools) have enabled the monitoring of cultural heritage areas that are threatened by natural disasters (e.g., climate change, earthquakes) and unplanned human activities (e.g., land use changes) [26][27][28].In the same context, the first launched satellite that was used in the field of archaeology, known as the Corona satellite, was produced for espionage and was launched in the 1960s.This satellite became a very important archive for archaeologists due to its high-resolution data and old landscape features [29].
Later, civilian satellites started with the launching of the family of Landsat (e.g., Land-sat-1 on 23 July 1972, Landsat-2 on 22 January 1975, Landsat-3 on 5 March 1978, and Landsat-4 on 16th July 1982), which continuously supported archaeologists with the required information with a spatial resolution data between 60 and 30 m [30].Based on the need for clear data, multi-series satellite systems, including the French SPOT HRV, Orbview (GeoEye), Quickbird, and Sentinel-2, which are characterized by their high spatial resolution, were launched.In addition, more recent hyperspectral sensors and radar satellites (e.g., Cosmo-Skymed, Sentinel-1) were installed and enabled workers in the field of archaeology to discover several new archaeological sites [31,32].
Ten years ago, the Egyptian Antiquities Authority (EAA) declared that the Pyramids Plateau was threatened by shallow groundwater levels close to the Great Pyramids and the Sphinx in Giza, Egypt.Thus, Egypt and USAID worked together to lower the groundwater to safe levels to protect the archaeological site from further damage [33].Unfortunately, recent studies revealed high groundwater levels close to the archaeological area [34], which has led to the need for additional extensive and serious work to provide solutions to reduce the groundwater level close to this heritage site.
Additionally, multiple geophysical measurements have been carried out in the region of the Great Pyramids in Giza and have proven that the area suffers from the nearsubsurface groundwater level, which strongly affects some archaeological sites in the studied area [35] (Figure 1a).The Giza archive documents show not only a 4500-year-old urban history but also a rapid stream of environmental events [36].Along with long and se-Land 2024, 13, 471 3 of 25 quential environmental events, the phenomenon of urban growth had undesirable impacts on the entire Cairo region, particularly in the heritage area of the Great Giza pyramids close to Nazlet El-Samman city [37].Indeed, geophysics measurements of the Giza Pyramids have indicated that the elevation of the groundwater table ranges mostly between 14.5 m and 17 m in the Khafre Causeway, Sphinx, and Sphinx Temple, central field of Mastaba, Valley Temple of Khafre, and Nazlet El-Samman areas [38].In Nazlet El-Samman village, geophysical surveys revealed that the elevation of the groundwater was approximately 17 m.This elevation mostly has a link with the Nile alluvium aquifer, which might raise the water table level below the Sphinx area, causing a large risk to the heritage site.The groundwater values of the measured sites revealed that the elevation was between 15 m and 15.5 m.The base-level elevation of the Sphinx is 20 m lower than that of the suburban area of Nazlet El-Samman, which might indicate a recharge of the aquifer below the Sphinx and increased capillary water rise.In the Sphinx area, the groundwater table based on different techniques has an elevation of approximately 15.5 m to 15.8 m.In contrast, the Valley Temple of Khafre and the central field of Mastaba measured a groundwater elevation of approximately 15 m.On the other hand, the Queen Khentkawes tomb has a groundwater elevation of approximately 14.8-15.2m.Similarly, the Menkaure Valley Temple had an average groundwater elevation of approximately 16 m (Figure 1b).According to the landscape, environment, and geo-environmental status of the study area, great attention should be given to the built-up spaces to the west of the Great Pyramids of Giza, which undoubtedly affected the groundwater table elevation of the area.Thus, a dewatering system should be implemented to lower the groundwater elevation to avoid aggressive hazards at archaeological sites [38].
Land 2024, 13, x FOR PEER REVIEW 3 of 24 4500-year-old urban history but also a rapid stream of environmental events [36].Along with long and sequential environmental events, the phenomenon of urban growth had undesirable impacts on the entire Cairo region, particularly in the heritage area of the Great Giza pyramids close to Nazlet El-Samman city [37].Indeed, geophysics measurements of the Giza Pyramids have indicated that the elevation of the groundwater table ranges mostly between 14.Similarly, the Menkaure Valley Temple had an average groundwater elevation of approximately 16 m (Figure 1b).According to the landscape, environment, and geo-environmental status of the study area, great attention should be given to the built-up spaces to the west of the Great Pyramids of Giza, which undoubtedly affected the groundwater table elevation of the area.Thus, a dewatering system should be implemented to lower the groundwater elevation to avoid aggressive hazards at archaeological sites [38].
(a) (b) Remotely sensed land and groundwater indicators are more capable of providing useful data and information than traditional techniques.Thus, integrated remote sensing and GIS applications are widely utilized in groundwater mapping and have become the main focus of numerous geo-archaeological research, beginning with aerial pictures [40].Therefore, this paper aims to assess the current status of the Giza heritage area based on the integration between remote sensing and GIS techniques.In addition, the study focuses on identifying and mapping the areas that are influenced by uncontrolled urban expansion and land use change, along with land subsidence.To achieve the study goals and quantify urban sprawl over the years, historical multi-temporal satellite data made Remotely sensed land and groundwater indicators are more capable of providing useful data and information than traditional techniques.Thus, integrated remote sensing and GIS applications are widely utilized in groundwater mapping and have become the main focus of numerous geo-archaeological research, beginning with aerial pictures [40].Therefore, this paper aims to assess the current status of the Giza heritage area based on the integration between remote sensing and GIS techniques.In addition, the study focuses on identifying and mapping the areas that are influenced by uncontrolled urban expansion and land use change, along with land subsidence.To achieve the study goals and quantify urban sprawl over the years, historical multi-temporal satellite data made up of Corona and Landsat have been investigated along with the new acquisitions provided by Orbview (GeoEye), Sentinel-1, and Sentinel-2 data.Results provided quantitative information on changes in the uncontrolled urban expansion that occurred around Egyptian monumental sites.On the other hand, the changes in the topography of the study area (land subsidence) were calculated based on the radar Sentinel-1 and SRTM data to show the effects of land use change on the paleo-landscape setting.Finally, based on this remote sensing and GIS-based analysis of multi-temporal and multi-sensor data analysis, a mitigation strategy is proposed for supporting the preservation of the archaeological area.

The Study Area
This research focused on the Giza complex situated at 31 • 7 46.369E and 29 • 58 30.734N, which includes the Great Pyramids, funerary, and temples.This area is situated on a plateau approximately 66 m above sea level [41].The study area includes the Great Pyramid (the pyramid of "Khufu"), which was built in the mid-20th century BC for the fourth-dynasty pharaoh; the pyramid of his son, "Khafre" (the middle one in size), is slightly smaller than that of Khufu; and the pyramid of "Menkaure" (his grandson), which is the third smallest in size and forms a complex beside the Sphinx [42].This area is famed for the presence of various temples and tombs dated primarily to the Old Kingdom and later Egyptian Kingdoms, particularly the New until the Late Kingdoms (Pharaoh's age), in ancient Egyptian history [43,44] (Figure 2a,b).
Land 2024, 13, x FOR PEER REVIEW 4 of 24 up of Corona and Landsat have been investigated along with the new acquisitions provided by Orbview (GeoEye), Sentinel-1, and Sentinel-2 data.Results provided quantitative information on changes in the uncontrolled urban expansion that occurred around Egyptian monumental sites.On the other hand, the changes in the topography of the study area (land subsidence) were calculated based on the radar Sentinel-1 and SRTM data to show the effects of land use change on the paleo-landscape setting.Finally, based on this remote sensing and GIS-based analysis of multi-temporal and multi-sensor data analysis, a mitigation strategy is proposed for supporting the preservation of the archaeological area.

The Study Area
This research focused on the Giza complex situated at 31°7′46.369″E and 29°58′30.734″N, which includes the Great Pyramids, funerary, and temples.This area is situated on a plateau approximately 66 m above sea level [41].The study area includes the Great Pyramid (the pyramid of "Khufu"), which was built in the mid-20th century BC for the fourth-dynasty pharaoh; the pyramid of his son, "Khafre" (the middle one in size), is slightly smaller than that of Khufu; and the pyramid of "Menkaure" (his grandson), which is the third smallest in size and forms a complex beside the Sphinx [42].This area is famed for the presence of various temples and tombs dated primarily to the Old Kingdom and later Egyptian Kingdoms, particularly the New until the Late Kingdoms (Pharaoh s age), in ancient Egyptian history [43,44] (Figure 2a,b).

Material
Multiple satellite sensors were utilized in this research, including the optical satellite imagery Corona (KH-4A) captured on 25 January 1965, Landsat5 (TM) captured on 07 May 2009, and the radar SRTM (1_arc-second) satellite image captured on 22 September 2014.These free-charge data were downloaded from the U.S. Geological Survey (USGS) website [45] and were used to detect and estimate the changes in the built-up areas close to the study area.In addition, it enables the creation of the SWAT model derived from SRTM data (ArcSWAT 2012.10.6.24)obtained from [46].Additionally, the optical satellite Sentinel-2A image obtained from ESA [47] was captured on 25 December 2016, along

Material
Multiple satellite sensors were utilized in this research, including the optical satellite imagery Corona (KH-4A) captured on 25 January 1965, Landsat5 (TM) captured on 7 May 2009, and the radar SRTM (1_arc-second) satellite image captured on 22 September 2014.These free-charge data were downloaded from the U.S. Geological Survey (USGS) website [45] and were used to detect and estimate the changes in the built-up areas close to the study area.In addition, it enables the creation of the SWAT model derived from SRTM data (ArcSWAT 2012.10.6.24)obtained from [46].Additionally, the optical satellite Sentinel-2A image obtained from ESA [47] was captured on 25 December 2016, along with the radar Sentinel-1A (SLC) image acquired in 2016 and 2019 (Table 1).These data were used for calculating the changes in the built-up areas around the study area in addition to the effects of these changes on the study area topography.The data were analyzed, and the required information was obtained using SNAP 9.0, ArcMap 10.6.1, and Envi 5.3 software.

Methods
In this research, an integration between optical and radar data was implemented to detect the changes in the urban areas between 1965 and 2019, along with the geoenvironmental status of the region.The optical Corona, Landsat, Sentinel-2, and radar Sentinel-1 data were preprocessed, processed, and post-processed to measure the changes in the built-up area.Based on the vector built-up data extracted from the raster satellite imagery, some spatial statistics were processed to calculate the changes in the space values between 1965 and 2019 using ArcMap software 10.6.1.Furthermore, the radar SRTM data were used to evaluate the topographic setting of the study area.Thus, the land subsidence rate was calculated based on the changes in elevation from 6-30 July 2016.The following sections will explain the applied techniques in this study in detail.

Optical Data Processing
The optical Corona, Landsat-5, and Sentinel-2 satellite data were downloaded, preprocessed, and processed to extract the urban areas.The Corona data were georeferenced and geometrically corrected based on the coordinates of the study area (WGS_1984_UTM_Zone_36N).The Sentinel-2 image was resampled, the Landsat image was preprocessed, and all the data were clipped based on the geographical boundary of the study area.In addition, unsupervised classification was performed to accurately detect the urban layer, and finally, the data were converted from raster to polygon.Based on the vector built-up data extracted from the raster satellite imagery, the changes in the urban areas between 1965 and 2019 were calculated.

•
Sentinel-1 (SLC) for extracting urban areas; The Sentinel-1 radar image was analyzed to detect and calculate the area of builtup zones in 2019 based on a model arranged inside SNAP software (GraphBuilder tool) that includes eight principles steps as follows: Read, TOPSAR split, Apply Orbit File, Calibration, TOPSAR Deburst, Multilook, Terrain correction, and Write operators.The S-1 TOPS split is applied to create a new product with only a sub-swath covering the ROI.When the study area covers more than one sub-swath, the chosen sub-swaths must be merged by the S-1 TOPS Merge operator [48].Thus, to apply orbit information, the precise orbit files are automatically downloaded for Sentinel-1 products by SNAP and added to its metadata, such as the accurate position of the satellite during the acquisition of the data.After that, the satellite image was radiometrically calibrated, and the sub-swaths were debuted to remove the seamlines between the single bursts [49].The next step is to multilook the calibrated data.Multilooking is the process of determining which pixels are averaged in the range and/or azimuth direction in an SLC dataset.On the other hand, averaging also reduced some speckling effects and smoothed the image.In addition, because the study used a single image, a co-registration step was not required [50].After performing the multilook step, the terrain correction method is applied to remove the misleading effect of topography on the backscatter values to remove the geolocation errors [51].Finally, the processed data are written to save the output by the Write operator [52] in a tiff formatting to be readable in the ENVI and ArcMap software.

•
Topographic setting of the study area based on the SRTM data; The Soil and Water Assessment Tool (SWAT) is an open-source model installed in ArcMap software [53] that is used to classify and divide a watershed into subbasins according to the SRTM (digital elevation model), and the hydrologic parameters such as the area, slope, and length of the extracted subbasins are extracted.In SWAT, subbasins are subdivided into hydrologic response units (HRUs) consisting of homogeneous topographical, land use, and soil characteristics [54,55].In this study, SWAT Version 2012, published by the Agricultural Research Service and the Texas Agricultural Experiment Station, was used.The ArcSWAT ArcGIS extension evolved from AVSWAT2000, an ArcView extension developed for an earlier version of SWAT [56].In our case, the Nash-Sutcliffe coefficient (NSE) was used to validate the SWAT model (Equation ( 1)): where P_x represents the model outputs derived from other DEM data, P _0 represents the mean observed streamflow, P_0 is the observed streamflow, and n represents the number of measured streamflows [57].
• Sentinel-1 (SLC) for detecting land subsidence values; The Sentinel-1 (SLC) imagery dated 6 July, 30 July, and 15 December 2016 were considered for calculating the land subsidence in the region over the short-term (between 6 and 30 July 2016) and long-term (between 30 July and 15 December 2016).The input data contained two single-look complex (SLC) data covering the study area in zip file format and were analyzed using SNAP software to perform the interferometric step.Thus, the split and apply orbit steps were applied by reading the two chosen data products (6 and 30 July-30 July and 15 December), specifically the VV polarization, selecting a sub-swath with TOPSAR-Split, and choosing the study area.After that, the two images were co-registered into a stack.The image on the first date (6 July-30 July 2016) was selected as the master image, and the other image (30 July-15 December 2016) was selected as the slave image, whereas each ground target contributed to both images having the same (range, azimuth) pixel.Thereafter, based on a series of steps automatically occurring based on a graph prepared on the SNAP software, the TOPS co-registered.Thus, the TOPS Deburst operator was used to seamlessly join all bursts in a swath into a single image.Subsequently, to emphasize phase signatures related to deformation, topographic phase contributions are typically removed using a known DEM.Then, to be able to properly analyze the phase signatures in the interferogram, the signal-to-noise ratio will be increased by applying multilooking and phase filtering techniques.In addition, for multilooking, a phase filtering step using a state-of-the-art filtering approach is performed.After the phase filtering step, the interferometric phase is significantly improved, and the value of the land deformation is clearly visible [58][59][60].Finally, the output data were exported to GeoTIFF format to match the rate of land subsidence phenomena with other factors (e.g., urban sprawl and the SWAT model).

Spatial Statistical Analysis
The data obtained from the satellite was analyzed for spatial patterns using various techniques such as cluster outliers and standard distance, geographically weighted regression, K function multidistance spatial clustering, optimized hot spot analysis, and spatial autocorrelation by distance.These techniques were applied to the spaces occupied by built-up areas.

•
Clusters outliers and standard distance; The process mechanism of this spatial statistics method is that the cluster and outlier analysis classified the processed layer (feature) into various groups or anomalous values based on the proximity criterion (attributes) to five degrees of geographical classes.The types of classification are (i) the spots that have either high or (ii) low values in concordance with their surroundings, (iii) anomalous areas where a spot has a value that is very different from the neighbors in higher or (iv) lower value, or (v) no relations can be made.This kind of process is mathematically implemented in ArcMap as follows (Equations ( 2) and ( 3)): where x_i represents an attribute for a feature i, X represents the mean of the corresponding attribute, x_(i, j) represents the spatial weight between features i and j, and n represents the total number of features as in Equation (3) [61]: • Spatial autocorrelation statistic using geographically weighted regression; Geographically weighted regression (GWR) is a familiar tool that can easily produce detailed maps of spatial variations in relationships.In addition, geographically weighted regression is a tool that is used for exploring spatial heterogeneity [62].This method also has the benefit of being able to produce local measures of conditional spatial dependency [63].As in Equation ( 4), it seems intuitively appealing to base estimates a_ik on observations taken at sample points close to i.If some degree of smoothness of a_iks is assumed, then reasonable approximations may be made by considering the relationship between the observed variables in a region geographically close to i as follows: where a_ik is the value of the kth parameter at location i [64].GWR allows coefficients to vary continuously over the study area, and a set of coefficients can be estimated at any location-typically on a grid-so that a coefficient surface can be visualized and interrogated for relationship heterogeneity.In essence, GWR measures the inherent relationships around each regression point i, where each set of regression coefficients is estimated by weighted least squares.The mathematical expression for this case is explored in Equation ( 5 where X is the matrix of the independent variables with a column of 1 s for the intercept; y is the dependent variable vector; βi = (β i0 , . . .β im ) T is the vector of m + 1 local regression coefficients.W i is the diagonal matrix denoting the geographical weighting of each observed data for regression point i [65].
• Space-time Ripley's K Function for Spatiotemporal Point Pattern Analysis Space-time Ripley's K function is a statistical approach that computes space-time point events and estimates their second-order property based on point pair distance calculations [66].It takes both the number of points and distances between points into account, which allows for quantitatively evaluating how much the observed point pattern deviates from randomness at multiple spatiotemporal scales.The theoretical space-time Ripley's K function is calculated by dividing E, the expected number of points within spatial distance s and temporal distance t, by the point intensity λst: Equation ( 6) characterizes the pattern of spatiotemporal points [67], where a cylinder of base πs2 and height 2t is centered on each point to count the number of neighboring points falling within.Then, the total number of points n is divided by the volume of the irregular cylinder formed by the study area and study duration, resulting in λst.It is expected that (s, ) = 2πs2t, if the point distribution obeys complete spatiotemporal randomness (CSTR), K(s,t) > 2πs2t if the points fit clustering within spatial distance s and temporal distance t, and K(s, t) < 2πs2t for dispersed space-time patterns.Using Equation (7), the space-time K function is formulated as: where A denotes the area of the study region and D is the duration of the study period.The product of A and D results in the volume of the irregular cylinder that is formed by the study area (base) and study period (height).ωij and υij are spatial and temporal weighting functions that correct edge effects respectively.dij is the spatial distance between points i and j, and uij is the temporal distance between points i and j.
The statistics that allow the expression of the value of a variable taken in the neighborhood are defined by the product of a measurement allowing the establishment of the similarity, or dissimilarity, cij, between the values of y, and the indicator of the spatial connectivity between the variables, wij (Equation ( 8)): where ωij is the element located on line i and column j of the spatial weights matrix (W).The measurement of similarity or dissimilarity, cij, usually allows us to gauge the gap between the values taken in the neighborhood, yj and the value of the variable considered, yi.The statistics of the products are then normalized by various multiplicative constants.There are different measurements of similarity/dissimilarity between the values of variable y.The different measurements of detection are based on these differences in the specification of the measurements used [68].

• Optimized Hot Spot Analysis
The Getis-Ord Gi * statistic [69][70][71][72] is a local statistical model that quantifies the degree of spatial dependence on the pronounced scale distance derived above.The Gi* statistics were implemented to determine the clustering level of the PS/DS MP, specifying a single MP at location i and its neighbors j within a searching distance d.For each single MP at site i, the Gi * index can be estimated as: Here, n is the total number of PS/DS MPs, nij is the number of PS/DS MPs within the scale distance d, x is the surface displacement recorded by PS/DS, and x and s are the mean value and the standard deviation of all displacements, respectively, that are revealed by the whole PS/DS MP.Through these Getis-Ord Gi* statistics, clustered and random patterns can be separately identified.The outcomes of the Getis-Ord Gi* statistics for each PS/DS MP are the calculated values of the z score (standard deviation) and p value (independence probability), both of which are utilized to estimate the statistical significance of spatial autocorrelation.Positive and negative z scores indicate spatial clusters of PS/DS moving towards and away from the LOS, respectively.Figure 3 shows the data collected, processing steps, outputs, and the proposed solution of the study area problem.
the data collected, processing steps, outputs, and the proposed solution of the study area problem.

Detect the Built-Up Area Time-Series Changes
The total area of the study area is estimated to be 66,614 km 2 , and the urban area represents 5083 km 2 of this total area, which is equivalent to 7.63% of the total area of the study area in 1965.In 2009, the urban mass area reached 21,798 km 2 , accounting for 32.72% of the total area of the study, with an annual increase of 0.379 km 2 /year.The urban mass area in 2016 was 31,087 km 2 , approximately half of the total area, with an annual increase of 1.327 km 2 /year between 2009 and 2016.In 2019, the urban mass area reached 39,782 km 2 , which is nearly two-thirds of the total area, with a high annual growth rate of approximately 2.898 km 2 /year.The annual growth rate between 1965 and 2019 was approximately 0.642 km 2 /year.The urban areas expanded approximately eight times in nearly 54 years and doubled twice in the last 10 years between 2009 and 2019.This can be attributed to significant population growth and security instability following the January 2011 revolution, which marked a turning point in continuous urban sprawl.This led to the establishment of new urban settlements that were not previously present in this area, as shown in Figure 4 and Table 2.

Detect the Built-Up Area Time-Series Changes
The total area of the study area is estimated to be 66,614 km 2 , and the urban area represents 5083 km 2 of this total area, which is equivalent to 7.63% of the total area of the study area in 1965.In 2009, the urban mass area reached 21,798 km 2 , accounting for 32.72% of the total area of the study, with an annual increase of 0.379 km 2 /year.The urban mass area in 2016 was 31,087 km 2 , approximately half of the total area, with an annual increase of 1.327 km 2 /year between 2009 and 2016.In 2019, the urban mass area reached 39,782 km 2 , which is nearly two-thirds of the total area, with a high annual growth rate of approximately 2.898 km 2 /year.The annual growth rate between 1965 and 2019 was approximately 0.642 km 2 /year.The urban areas expanded approximately eight times in nearly 54 years and doubled twice in the last 10 years between 2009 and 2019.This can be attributed to significant population growth and security instability following the January 2011 revolution, which marked a turning point in continuous urban sprawl.This led to the establishment of new urban settlements that were not previously present in this area, as shown in Figure 4 and Table 2.

Detect the Changes by Clusters Outliers and Standard Distance
Based on Figure 5, it can be observed that the standard distances in shapes a, b, c, and d vary and are inconsistent.The center of the circle changes corresponding to the increase in urban mass area in the study area between 1965 and 2019.Shapes a and b are more concentrated and clustered as their circle radius decreases, while shapes c and d are more dispersed, indicating the spread of urban masses as their circle radius increases.This observation is further supported by the locations of the three pyramids.At the beginning of the period (Figure 5a), they were situated far from the circle representing the standard distance, while at the end of the period, they were closer to the central point of the figure representing the standard distance in 2019 (Figure 5d).
For the cluster outliers, it can be observed that the urban masses in the years 1965-2009 exhibit random distributions, as well as dispersion and divergence.The areas with a Z score less than 1, especially those ranging between (1: −1), whether in the western, northern, or southern regions, are shown in Figure 5a,b.The distances between urban masses converge and cluster in the center of the area in small spaces, while the remaining urban masses spread over a large area on the outskirts of the region.In contrast, in 2016 and 2019 (Figure 5c,d), an increase in the areas with a Z score greater than one, especially those with a Z score exceeding two, was observed.This indicates that the dominant urban pattern is clustering and convergence, especially in the eastern and western regions of the study area.This is due to the increase in the urban mass area in the study area.

Detect the Changes by Clusters Outliers and Standard Distance
Based on Figure 5, it can be observed that the standard distances in shapes a, b, c, and d vary and are inconsistent.The center of the circle changes corresponding to the increase in urban mass area in the study area between 1965 and 2019.Shapes a and b are more concentrated and clustered as their circle radius decreases, while shapes c and d are more dispersed, indicating the spread of urban masses as their circle radius increases.This observation is further supported by the locations of the three pyramids.At the beginning of the period (Figure 5a), they were situated far from the circle representing the standard distance, while at the end of the period, they were closer to the central point of the figure representing the standard distance in 2019 (Figure 5d).
For the cluster outliers, it can be observed that the urban masses in the years 1965-2009 exhibit random distributions, as well as dispersion and divergence.The areas with a Z score less than 1, especially those ranging between (1: −1), whether in the western, northern, or southern regions, are shown in Figure 5a,b.The distances between urban masses converge and cluster in the center of the area in small spaces, while the remaining urban masses spread over a large area on the outskirts of the region.In contrast, in 2016 and 2019 (Figure 5c,d), an increase in the areas with a Z score greater than one, especially those with a Z score exceeding two, was observed.This indicates that the dominant urban pattern is clustering and convergence, especially in the eastern and western regions of the study area.This is due to the increase in the urban mass area in the study area.

Detect the Changes by Geographically Weighted Regression
Based on Figure 6, the weighted geographic regression analysis indicated that the study area exhibited spatial heterogeneity during the period from 1965 to 2019.The red areas on the map indicate higher rates of urban expansion compared to the four forms (a,

Detect the Changes by Geographically Weighted Regression
Based on Figure 6, the weighted geographic regression analysis indicated that the study area exhibited spatial heterogeneity during the period from 1965 to 2019.The red areas on the map indicate higher rates of urban expansion compared to the four forms (a, b, c, and d).These areas extended significantly in the inner parts of the study area during 1965, 2009, and 2016, and spread to the northern and southern parts in 2019.On the other hand, the blue areas on the map represent regions that experienced more urban growth than expected.Despite not being expected to occur, urban clusters spread significantly on the margins of the northern, eastern, southern, and western parts of the study area during the entire study period.
hand, the blue areas on the map represent regions that experienced more than expected.Despite not being expected to occur, urban clusters spread s the margins of the northern, eastern, southern, and western parts of the s ing the entire study period.

Multidistance Spatial Cluster (K Function) Time-Series Changes
The mathematical transformation of Ripley's K-function performed slightly different, and as a final result, the expected index of a random patt the input distance [64].To assess the pattern type, the observed L(d) compared with the expected value.In each distance increment, the actual compared against a random distribution, and if the observed L(d) value expected, a clustered pattern is dominant at that scale of analysis.In turn, smaller than expected, the distribution is considered to be dispersed.Mo index of clustering or dispersion to be statistically significant, it must smaller, respectively, than the high or low level of the confidence interv study, Figure 7a indicates that there is a clustering and agglomeration of

Multidistance Spatial Cluster (K Function) Time-Series Changes
The mathematical transformation of Ripley's K-function performed in ArcGIS is slightly different, and as a final result, the expected index of a random pattern is equal to the input distance [64].To assess the pattern type, the observed L(d) value must be compared with the expected value.In each distance increment, the actual distribution is compared against a random distribution, and if the observed L(d) value is greater than expected, a clustered pattern is dominant at that scale of analysis.In turn, if this value is smaller than expected, the distribution is considered to be dispersed.Moreover, for the index of clustering or dispersion to be statistically significant, it must be greater or smaller, respectively, than the high or low level of the confidence interval [65].In our study, Figure 7a indicates that there is a clustering and agglomeration of urban masses that are concentrated near distances of approximately 800 to 2000 m in which the red line (observed) deviates above the blue line (expected).This deviation reached a maximum distance between the two lines at 1400 m, where it reached 600 concerning L(d).Additionally, Figure 7b indicates that there is dispersion and divergence in the urban masses when the distance is less than 400 m, where the blue line deviates above the red line.This dispersion and divergence gradually decrease until there is a concentration of urban masses, especially at distances ranging between 1200 and 2400 m.

Spatial Autocorrelation Statistic Time-Series Changes
In 1965, the Global Moran s Index Summary by Distance for the Spatial Autocorrelation statistic for the urban layer indicated that the beginning distance was 829, the distance increment was 13.95, the first peak (distance, value) was 842.95 and 2.53, respectively, and the maximum peak (distance, value) was 940.63 and 2.56, respectively.The results revealed that Moran s index is greater than zero, which indicates that the distribution is regular and clustered but not to a large extent (Table 3 and Figure 8a).All the Z score values were limited to the range of 1.96-2.58,and the highest value was 2.56 at a distance of 940.63 m, while the lowest value was 2.23 when the distance was 884.81 m.Furthermore, the lowest distance between urban masses was nearly 829 m, while the largest distance was 954.58 m.This closeness between the distances indicates that this clustering is regular, which reflects the unfair distribution of urban masses in the region in 1965.

Distance (m)
Moran's Index z-Score Figure 7c indicates the dispersion and separation of urban masses when the distance is less than approximately 500 m, and there was an agglomeration of urban blocks, especially when the distance exceeded 1400 m. Figure 7d shows that dispersion and separation increased to reach distances less than 600 m, while there was urban clustering and agglomeration at distances exceeding 600 m.The deviation between the red and blue lines decreased, and they became closer to each other, as shown in Figure 7c,d, as a result of the increase in urban growth between 2016 and 2019.The urban masses moved closer together due to increasing urban growth.

Spatial Autocorrelation Statistic Time-Series Changes
In 1965, the Global Moran's Index Summary by Distance for the Spatial Autocorrelation statistic for the urban layer indicated that the beginning distance was 829, the distance increment was 13.95, the first peak (distance, value) was 842.95 and 2.53, respectively, and the maximum peak (distance, value) was 940.63 and 2.56, respectively.The results revealed that Moran's index is greater than zero, which indicates that the distribution is regular and clustered but not to a large extent (Table 3 and Figure 8a).All the Z score values were limited to the range of 1.96-2.58,and the highest value was 2.56 at a distance of 940.63 m, while the lowest value was 2.23 when the distance was 884.81 m.Furthermore, the lowest distance between urban masses was nearly 829 m, while the largest distance was 954.58 m.This closeness between the distances indicates that this clustering is regular, which reflects the unfair distribution of urban masses in the region in 1965.Additionally, for 2009, the Global Moran s Index Summary by Distance for the Spatial Autocorrelation statistic for the urban layer showed that the initial distance was 541.0, the distance increment was 77.8, the first peaks (distance, value) were 696.6 and 22.1, respectively, and the maximum peaks (distance, value) were 1007.8 and 25.77, respectively.The results indicated that the Moran coefficient increased compared to that in the previous period, and the z score also increased significantly, ranging between 19.06 when the distance was 541 m and 25.77 when the distance was 1007.82 m (Table 4 and Figure 8b).This demonstrates the agglomeration and clustering of urban masses due to the great urban expansion witnessed by the region between 1965 and 2009.This leads to all Z score Additionally, for 2009, the Global Moran's Index Summary by Distance for the Spatial Autocorrelation statistic for the urban layer showed that the initial distance was 541.0, the distance increment was 77.8, the first peaks (distance, value) were 696.6 and 22.1, respectively, and the maximum peaks (distance, value) were 1007.8 and 25.77, respectively.The results indicated that the Moran coefficient increased compared to that in the previous period, and the z score also increased significantly, ranging between 19.06 when the distance was 541 m and 25.77 when the distance was 1007.82 m (Table 4 and Figure 8b).This demonstrates the agglomeration and clustering of urban masses due to the great urban expansion witnessed by the region between 1965 and 2009.This leads to all Z score values increasing to more than 2.58.For 2016, the Global Moran's Index Summary by Distance for the Spatial Autocorrelation statistic for the urban layer proved that the beginning distance was 445, the distance increment was 86.22, and the first peak (distance, value) and max peak (distance, value) were recorded as none.The results revealed that the Moran coefficient increased compared to that in the previous period, as did the z score.The z score ranged between 15.94 when the distance was 445 m and 28.70 when the distance was 1220.98 m (Table 5 and Figure 8c).This indicated that urbanization became closer and more clustered during this period in 2016 than during the previous period in 2009.In 2019, the Global Moran's I Summary by Distance for the Spatial Autocorrelation statistic for the urban layer showed that the initial distance was 702, the distance increment was 91.07, the first peaks (distance, value) were 1066.28 and 37.89, respectively, and the maximum peaks (distance, value) were 1248.41 and 38.83, respectively.The results in Table 6 and Figure 8d show that the distance between urban masses increased due to the newly developed urban communities along the region's outskirts.Additionally, a noticeable increase in the z score was detected, ranging between 31.47 when the distance was 793.07 m and 38.83 when the distance was 1248.41 m.This stage is considered the urban explosion in the region, in which urbanization became more clustered, agglomerated, and closer in this period of 2019 compared to the previous periods of study.The study area was analyzed using Sentinel-1 (SLC) data that were collected over the short and long terms for showing the effects of urban expansion on the rise in groundwater levels in the area around the archaeological area of Giza Necropolis.For instance, Bokhari et al. (2023) analyzed land subsidence based on long-term data spanning from July 2017 to July 2019 [73].Similarly, Ghazifard et al. (2017) evaluated subsidence rates during two extended periods, from March to December 2005 and from July 2011 to January 2012.Their findings indicated that the primary causes of ground subsidence were declining groundwater levels and changes in sediment composition and thickness [74].Additionally, Tiwari et al. (2017) estimated subsidence over an urban area using a time series of Sentinel-1 data from 2014 to 2017, with an image acquired on 29 May 2016, serving as the master image, resulting in a maximum perpendicular baseline of 150 m [75].Bhattarai et al. (2017) detected land subsidence in the Kathmandu Valley over an extended period from 2007 to 2010 [76].Furthermore, Nalakurthi and Behera (2022) examined the rate of land subsidence in Mumbai and investigated the impact of groundwater exploitation using Sentinel-1 SAR imagery collected from October 2014 to May 2019.Mumbai, being one of the most densely populated cities in the world with significant urban sprawl and water scarcity issues, experiences groundwater exploitation for domestic and industrial purposes, resulting in land subsidence and affecting coastal communities due to the influence of sea levels [77].
In this study, the two chosen dates were between 6 and 30 July 2016 and in the period between 30 July and 15 December 2016.The data showed that land subsidence occurred in this area, as illustrated in Figure 9a,b.The changes were mostly observed in the southern region near Nazlet El-Samman, which is close to the pyramids.On the other hand, the northern region experienced little change.The subsidence in the southern region reached about −0.0138 m over the period between 30 July and 15 December 2016.It is believed that groundwater may have caused the land subsidence, leading to the gradual sinking or settling of the land surface.When water seeps into the ground and saturates the soil layers, the soil can lose its structural integrity and contract.Consequently, structures such as the Great Pyramids and the Sphinx may have started sinking, which could eventually lead to structural damage.
The historical geophysical characteristics showed that there was a waterway in this region, with a groundwater level ranging between 14.2-16.3m, and higher values were observed close to the Pyramid of Khufu, which is approximately 17.4-19.4m [38].Moreover, urban expansion increased in recent years, particularly in 2019, which put pressure on the area and increased land subsidence.The threat facing the Great Pyramids and Sphinx region arose from the impact of groundwater leakage from nearby irrigation canals and condensed urban masses [35].
The historical geophysical characteristics showed that there was a waterway in this region, with a groundwater level ranging between 14.2-16.3m, and higher values were observed close to the Pyramid of Khufu, which is approximately 17.4-19.4m [38].Moreover, urban expansion increased in recent years, particularly in 2019, which put pressure on the area and increased land subsidence.The threat facing the Great Pyramids and Sphinx region arose from the impact of groundwater leakage from nearby irrigation canals and condensed urban masses [35].Overall, the results show that urban expansion greatly influenced the study area, particularly in 2019.Urban masses increased and clustered in the region, which was proven by different indices, such as geographically weighted regression, the Multidistance spatial cluster (K function), and Moran s index.The study revealed that newly developed urban communities led to an increase in the spaces between urban masses along the region s outskirts in 2019.This new development in the region resulted in urban explosions becoming more clustered and condensed in 2019 compared to previous years.This unmanaged urban expansion resulted in increased water leakage to groundwater in neighboring areas.The great pyramids of Giza have been threatened by this rise in groundwater level resulting from water infiltration from nearby urban masses.Two main aquifers are located near the Sphinx, with water levels of 1.5-4 m and 4-7 m below the surface [34].These aquifers are recharged by water leakage from urban areas, leading to severe hazards in the area.Pyramid construction materials and structures might be damaged as a result of this water infiltration.

Recommendation
Monitoring and regulating groundwater levels, implementing efficient drainage systems, and managing leaks from urban growth and irrigation canals are essential to reduce their impact on the Great Pyramids and Sphinx area.To reduce the impact of groundwater on the Great Pyramids and the Sphinx area, the stabilization of structures and moisture management should also be given top priority in preservation and restoration activities.Protecting cultural heritage, such as archaeological and cultural sites, is undoubtedly considered a critical issue and an urgent need [78][79][80][81].In addition, the Overall, the results show that urban expansion greatly influenced the study area, particularly in 2019.Urban masses increased and clustered in the region, which was proven by different indices, such as geographically weighted regression, the Multidistance spatial cluster (K function), and Moran's index.The study revealed that newly developed urban communities led to an increase in the spaces between urban masses along the region's outskirts in 2019.This new development in the region resulted in urban explosions becoming more clustered and condensed in 2019 compared to previous years.This unmanaged urban expansion resulted in increased water leakage to groundwater in neighboring areas.The great pyramids of Giza have been threatened by this rise in groundwater level resulting from water infiltration from nearby urban masses.Two main aquifers are located near the Sphinx, with water levels of 1.5-4 m and 4-7 m below the surface [34].These aquifers are recharged by water leakage from urban areas, leading to severe hazards in the area.Pyramid construction materials and structures might be damaged as a result of this water infiltration.

Recommendation
Monitoring and regulating groundwater levels, implementing efficient drainage systems, and managing leaks from urban growth and irrigation canals are essential to reduce their impact on the Great Pyramids and Sphinx area.To reduce the impact of groundwater on the Great Pyramids and the Sphinx area, the stabilization of structures and moisture management should also be given top priority in preservation and restoration activities.Protecting cultural heritage, such as archaeological and cultural sites, is undoubtedly considered a critical issue and an urgent need [78][79][80][81].In addition, the preservation of cultural heritage should be planned as a means of promoting conservation, documentation, and even restoration based on a specific strategy [82][83][84][85].Therefore, many efforts have been made to address the issues of urban expansion and groundwater leakage and their massive impacts on the Great Pyramids and the Sphinx region, in addition to mitigating their impact on the region.Based on the topography, built-up areas, and subsidence value properties, this research suggested a scenario for minimizing the groundwater level close to the Giza Pyramids area.Figure 10a-d illustrates the factors that influence groundwater leakage in the study area, which includes the topography, urban mass distribution, and land subsidence level.Figure 10e shows the suitable locations of the proposed trenches that could help reduce the impact of the groundwater level on the Great Pyramids and the Sphinx.

Conclusions
Changes in land use and land cover can have significant impacts on archaeological areas such as the Giza Necropolis in Egypt.Detecting the impacts of land-use change on this archaeological area is essential for its preservation and management.Therefore, this study aimed to track the extent of land-use changes over time for 1965, 2009, 2016, and 2019 based on remote sensing and GIS integration.In addition, the relationships between groundwater leakage and land subsidence in the region were studied.Multiple satellite sensors were utilized in the current study such as Landsat, Corona, as well Sentinel 1 and 2. These satellite data were analyzed based on different techniques, including cluster outliers, the Moran index, and spatial autocorrelation, to examine the changes in urban masses between the 1965 and 2019 period.This study revealed that urban expansion and infrastructure development can, directly and indirectly, impact the Giza Necropolis.The land use changed in the study area due to human interactions, in which the urban masses increased from 5.083 km 2 in 1965 to 39.782 in 2019 with an annual growth rate of 0.642 Km 2 /year.Urbanization led to encroachment, increased population density, construction, and the establishment of modern settlements near the archaeological site.These changes resulted in physical damage, disturbance of underground structures, increased pollution, and alteration of the surrounding landscape.Understanding the linkage between land-use changes and archaeological areas enables the identification of areas at greater risk from urban sprawl and the identification of areas that require immediate attention for preservation efforts.This study proposed new locations for trenches that can minimize the effect of groundwater discharge and reduce the impact on the Giza Pyramids and Sphinx.

Conclusions
Changes in land use and land cover can have significant impacts on archaeological areas such as the Giza Necropolis in Egypt.Detecting the impacts of land-use change on this archaeological area is essential for its preservation and management.Therefore, this study aimed to track the extent of land-use changes over time for 1965, 2009, 2016, and 2019 based on remote sensing and GIS integration.In addition, the relationships between groundwater leakage and land subsidence in the region were studied.Multiple satellite sensors were utilized in the current study such as Landsat, Corona, as well Sentinel 1 and 2. These satellite data were analyzed based on different techniques, including cluster outliers, the Moran index, and spatial autocorrelation, to examine the changes in urban masses between the 1965 and 2019 period.This study revealed that urban expansion and infrastructure development can, directly and indirectly, impact the Giza Necropolis.The land use changed in the study area due to human interactions, in which the urban masses increased from 5.083 km 2 in 1965 to 39.782 in 2019 with an annual growth rate of 0.642 Km 2 /year.Urbanization led to encroachment, increased population density, construction, and the establishment of modern settlements near the archaeological site.These changes resulted in physical damage, disturbance of underground structures, increased pollution, and alteration of the surrounding landscape.Understanding the linkage between land-use changes and archaeological areas enables the identification of areas at greater risk from urban sprawl and the identification of areas that require immediate attention for preservation efforts.This study proposed new locations for trenches that can minimize the effect of groundwater discharge and reduce the impact on the Giza Pyramids and Sphinx.
5 m and 17 m in the Khafre Causeway, Sphinx, and Sphinx Temple, central field of Mastaba, Valley Temple of Khafre, and Nazlet El-Samman areas [38].In Nazlet El-Samman village, geophysical surveys revealed that the elevation of the groundwater was approximately 17 m.This elevation mostly has a link with the Nile alluvium aquifer, which might raise the water table level below the Sphinx area, causing a large risk to the heritage site.The groundwater values of the measured sites revealed that the elevation was between 15 m and 15.5 m.The base-level elevation of the Sphinx is 20 m lower than that of the suburban area of Nazlet El-Samman, which might indicate a recharge of the aquifer below the Sphinx and increased capillary water rise.In the Sphinx area, the groundwater table based on different techniques has an elevation of approximately 15.5 m to 15.8 m.In contrast, the Valley Temple of Khafre and the central field of Mastaba measured a groundwater elevation of approximately 15 m.On the other hand, the Queen Khentkawes tomb has a groundwater elevation of approximately 14.8-15.2m.

Figure 1 .
Figure 1.Groundwater level close to study area: (a) partial deterioration in low part of Sphinx body due to near-subsurface groundwater [39]; (b) groundwater elevation values in area of Great Pyramids with particular focus on area of the Sphinx (based on [38]).

Figure 1 .
Figure 1.Groundwater level close to study area: (a) partial deterioration in low part of Sphinx body due to near-subsurface groundwater [39]; (b) groundwater elevation values in area of Great Pyramids with particular focus on area of the Sphinx (based on [38]).

Figure 2 .
Figure 2. It shows the study area: (a) location map of Egypt highlighting area of interest obtained from Google Earth and Sentinel-2 imagery with the Great Giza Pyramids and Sphinx; (b) hillshade of study area including contour lines and elevation points obtained from SRTM 1 Arc-Second Global.

Figure 2 .
Figure 2. It shows the study area: (a) location map of Egypt highlighting area of interest obtained from Google Earth and Sentinel-2 imagery with the Great Giza Pyramids and Sphinx; (b) hillshade of study area including contour lines and elevation points obtained from SRTM 1 Arc-Second Global.

Figure 3 .
Figure 3. Schematic diagram of methodology adopted in this research.

Figure 3 .
Figure 3. Schematic diagram of methodology adopted in this research.

Figure 4 .
Figure 4. Variations in built-up area in study area for period of 1965-2019.

Figure 4 .
Figure 4. Variations in built-up area in study area for period of 1965-2019.

Figure 9 .
Figure 9. Changes in land subsidence values of study area over short and long terms: (a) between 6 and 30 July 2016; (b) between 30 July and 15 December 2016.

Figure 9 .
Figure 9. Changes in land subsidence values of study area over short and long terms: (a) between 6 and 30 July 2016; (b) between 30 July and 15 December 2016.

Figure 10 .
Figure 10.Recommended model for lowering groundwater level close to the Giza Pyramids area based on: (a) topographic setting; (b) built-up areas; (c,d) subsidence values properties; (e) including recommended locations for trenches, pipeline, and water recycling station.

Author Contributions:
The authors declare that the conceptualization of this scientific research has been performed by A.E. and R.L.; the methodology has been performed by A.E.; the software has been performed by A.E.; the validation has been performed by A.E. and R.L.; the formal analysis has been performed by A.E.; the investigation has been performed by A.E., N.Z., W.M. and E.H.; the resources have been performed by A.E., N.Z., W.M. and E.H.; the data curation has been performed by N.Z., W.M. and E.H.; the writing-original draft preparation has been performed by A.E., N.Z., W.M. and E.H.; and the writing-review and editing have been performed by A.E. and R.L.All authors have read and agreed to the published version of the manuscript.Funding: This research received no external funding.

Figure 10 .
Figure 10.Recommended model for lowering groundwater level close to the Giza Pyramids area based on: (a) topographic setting; (b) built-up areas; (c,d) subsidence values properties; (e) including recommended locations for trenches, pipeline, and water recycling station.

Table 1 .
Data utilized and satellite imagery specifications.

Table 3 .
Moran s Index summary by Distance for Spatial Autocorrelation statistic for urban layer in 1965.

Table 3 .
Moran's Index summary by Distance for Spatial Autocorrelation statistic for urban layer in 1965.

Table 4 .
Moran's Index summary by Distance for Spatial Autocorrelation statistic for urban layer in 2009.

Table 5 .
Moran's Index summary by Distance for Spatial Autocorrelation statistic for urban layer in 2016.

Table 6 .
Moran's Index summary by Distance for Spatial Autocorrelation statistic for urban layer in 2019.