Mapping buried paleogeographical features of the Nile Delta (Egypt) using the Landsat archive

The contribution highlights the use of Landsat spectral-temporal metrics (STMs) for the detection of surface anomalies that are potentially related to buried near-surface paleogeomorphological deposits in the Nile Delta (Egypt), in particular for a buried river branch close to Buto. The processing was completed in the Google Earth Engine (GEE) for the entire Nile Delta and for selected seasons of the year (summer/winter) using Landsat data from 1985 to 2019. We derived the STMs of the tasseled cap transformation (TC), the Normalized Difference Wetness Index (NDWI), and the Normalized Difference Vegetation Index (NDVI). These features were compared to historical topographic maps of the Survey of Egypt, CORONA imagery, the digital elevation model of the TanDEM-X mission, and modern high-resolution satellite imagery. The results suggest that the extent of channels is best revealed when differencing the median NDWI between summer (July/August) and winter (January/February) seasons (ΔNDWI). The observed difference is likely due to lower soil/plant moisture during summer, which is potentially caused by coarser-grained deposits and the morphology of the former levee. Similar anomalies were found in the immediate surroundings of several Pleistocene sand hills (“geziras”) and settlement mounds (“tells”) of the eastern delta, which allowed some mapping of the potential near-surface continuation. Such anomalies were not observed for the surroundings of tells of the western Nile Delta. Additional linear and meandering ΔNDWI anomalies were found in the eastern Nile Delta in the immediate surroundings of the ancient site of Bubastis (Tell Basta), as well as several kilometers north of Zagazig. These anomalies might indicate former courses of Nile river branches. However, the ΔNDWI does not provide an unambiguous delineation.

1 Introduction

The reconstruction of the paleo-topography and paleo-hydrography of the Nile Delta plays a central role in landscape-focused archeological investigations. Historical settlement activity was strongly linked to and influenced by the presence and dynamics of the ancient river network of the delta. The great importance of the river network for early settlements is underlined by the fact that larger cities were only found in the immediate vicinity of larger Nile branches, which were of outstanding importance for traffic and trade and met the basic need of water for agriculture and food security (Bietak, 1975). The earliest textual sources that allow, to a certain degree, the reconstruction of the branches of the Nile mainly date from the 5th century BCE to the 4th century CE and come from various Greek and Roman authors (including Herodotus, Diodorus, Strabo, and Ptolemy). These sources named and described the estuaries of the Nile and the landscape/riverscape of the delta in more or less detail; however, they generally do not allow clear localization in today's topographic context. Most records indicate seven main branches of the Nile, which were either named after the cities at their mouths or after an important city located on the respective arm (Bietak, 1975; Ginau et al., 2019).

Today's situation differs greatly, as only two main Nile branches still exist: the western Rosetta and the eastern Damietta arms. The fluvial landscape of the Nile Delta has, therefore, changed tremendously, indicating that the delta was a highly dynamic environment. The same is true for the time before the classical historiographers. In the time of ancient Egyptian culture, settlement patterns and overall economic, cultural, and religious processes were greatly impacted by changes to the hydrographic system of the alluvial plain (Butzer, 1976). Especially in recent decades, the importance of understanding the paleo-topographies and paleo-hydrographies of the Nile Delta at different times has been recognized, and several research projects have explored the topic. Geophysical and geoarcheological investigations have been carried out in different regions of the delta at various scales. Most of these surveys relied principally on the sedimentological analysis of core drillings (e.g., Sewuster and van Wesemael, 1987; Andres and Wunderlich, 1991; Wunderlich 1988, 1989; Stanley and Warne, 1993a, b; Stanley et al., 1996; Flaux et al., 2012; Marriner et al., 2012; Ginau et al., 2019) and/or geophysical measurements (El-Gamili et al., 1994; El Gamili et al., 2001; El-Mahmoudi and Gabr, 2009; Pennington and Thomas, 2016).

Additionally, remote sensing data, as well as topographic and historical maps, have been widely used as they allow area-wide investigation and sometimes even provide records on the past topographic situation that might not be visible today. Bietak (1975), for instance, reconstructed the course of some of the Nile branches of the delta based on the analysis of topographic maps and the location of ancient settlements. He mapped linear structures, which presumably represent natural embankments of former rivers (levees). Since the end of the 1980s, satellite images (e.g., Landsat, SPOT, Corona, RapidEye) have been increasingly used to detect old river courses in the Nile Delta (e.g., Wunderlich, 1989; Marcolongo, 1992; Moshier and El-Kalani, 2008; Wilson and Grigoropoulos, 2009; Trampier et al., 2013; Ginau et al., 2017). Along with this, digital elevation models (DEMs) from the Shuttle Radar Topography Mission (SRTM; Stanley and Jorstad, 2006), as well as high-resolution data acquired by the TanDEM-X mission, have been used (Ginau et al., 2019). Lately, Elfadaly et al. (2020) showed the use of various remote sensing sources (including Landsat) to identify potential former settlement areas in the northern delta.

The above-listed studies all use geospatial data to indicate features of the modern land surface (e.g., location and orientation of field boundaries, anomalies of land cover, discontinuities of surface morphology) that may be related to the paleogeomorphological setting. Although these datasets do not allow chronological information to be derived or remains of the paleo-landscape to be detected without ambiguities, they are very valuable to support geophysical research on the ground as they help to narrow down the survey area, or indicate new promising locations for fieldwork. However, this requires buried features to have a distinct surficial expression compared to their surroundings so they can be detected with remotely sensed imagery.

In many of the studies listed above, results on the location of former landform features were obtained using single images or a few remote sensing datasets. With the opening of the Landsat Archive in 2008, an extensive time series of multispectral satellite data has become available to the public, open and free of charge. This archive compiles all available Landsat data, beginning with the first acquisitions made by Landsat 1 in the early 1970s. However, analysis of the entire Landsat archive for a region of interest is difficult using an individual processing strategy and standard infrastructure. This is due to high storage and processing requirements that – considering the huge amount of data – must be undertaken automatically. New analysis methods can provide a solution here, e.g., cloud-based infrastructure and automated processing chains. Both are offered by the Google Earth Engine (GEE) (Gorelick et al., 2017), a free cloud-based service that grants access to the entire Landsat archive and offers comprehensive algorithms for data processing and feature derivation. This enables processing of the Landsat Archive for the entire Nile Delta, and therefore offers a new and still unexploited pool of remote sensing data that can indicate features and anomalies related to the paleogeomorphological setting. The high spatial resolution and the long temporal baseline offered by the archive make it possible to visualize and investigate even small changes and subtle differences in reflection properties

This contribution presents results on the analysis of the Landsat Archive for the entire Nile Delta between 1985 and 2019 in order to detect buried paleogeographical features. The focus of the research is on the detection of anomalies that might indicate hints of a former, potentially ancient, geomorphological setting. Special emphasis is on the identification of (partially) buried geziras and tells and on potential courses of former or abandoned Nile branches. Anomalies detected in the Landsat time series are compared to the TanDEM-X DEM, historical topographic maps of the Survey of Egypt, and to satellite imagery of the CORONA mission and recent high-resolution imagery.

2 Study area

The Nile Delta (Fig. 1), an alluvial plain in the north of Egypt, today covers an area of about 24 000 km² and is densely populated and intensively used for agriculture (Pennington et al., 2017; Fig. 1). It is the youngest in a long series of deltaic formations that probably date back to the Miocene (Butzer, 1976; Said, 1981). The landscape dynamics of the delta were controlled by natural factors such as tectonics and climate and sea level fluctuations up to the Middle Holocene, but human influences have played an increasingly important role since the Late Holocene (Pennington et al., 2017).

Figure 1 Study area and data availability: (a) TanDEM-X digital elevation model (DEM) with a resolution of 13 m by 13 m per pixel and (b) Landsat-8 true-color composite (RGB) for median surface reflectance between 2015 and 2019 and number of cloud-free Landsat observations between 1985 and 2019 for the months of (c) July and August and (d) January and February. Insets in (a, b) are centered on the excavation site of Bubastis (south of the modern city of Zagazig) and show the full resolution of the datasets. Landsat imagery courtesy of the US Geological Survey. TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

In detail, the geological structure and history of the Nile Delta are quite complex and show regional differences (cf. Andres and Wunderlich, 1991; see also the delta-spanning summary of Pennington et al., 2017). In simplified terms, the geology can be described as follows. At the beginning of the Holocene, the delta was largely covered by sandy to sandy-gravel deposits of the Mit Ghamr and Geziracover formations. The Mit Ghamr formation consists of numerous smaller units of different genetic origin of Pleistocene age, deposited during an earlier, interwoven Nile regime (Said, 1981). Up to ca. 8000 cal BP, the Nile arms eroded or redeposited this material. On some topographic heights known as “turtlebacks” or “geziras”, an aeolian rearrangement of the fluvial sands forming the Geziracover formation often occurred (Wunderlich, 1989; Pennington et al., 2017). The formation of the Holocene alluvial delta plain began in the Middle Holocene. From about 8000 cal BP onwards, hydrological changes and increased sediment supply in the “African humid period” led to high accumulation rates in the Nile Delta and to the development of swampy wetlands.

The area was characterized by a widely ramified river network with extensive flood plains. The bluish-black, organic-rich, silty-clayey to clayey-silty deposits with intercalated peat horizons of this phase belong to the Bilqas 2 formation, which was deposited between ca. 8000 and 6000 cal BP in large parts of the delta. The relief of the delta was slightly hilly at this time, with topographic highs in many places, especially on the edges of the delta, extending several meters above the flood plain. The coastline was further inland than today, especially in the east (Goiran et al., 2005; Pennigton et al., 2017). Between 6000 and 5500 cal BP, relative sea level rise decreased and conditions became increasingly arid. The accumulation rate in the delta therefore decreased strongly and the river landscape changed noticeably. In the southern and central part of the delta, the sediments of the Bilqas 1 formation (brown-gray in color, and less rich in organic material) were already present at this time, deposited in a deltaic landscape with much wider, well-drained flood plains and individual river courses (Pennigton et al., 2017). By around 3500 cal BP, the sediments of the Bilqas 2 facies had almost completely disappeared and were spilled by sediments of the Bilqas 1 facies during the Late Holocene. The geziras, which were a typical landscape feature in the delta before 6000 cal BP, therefore gradually became rarer and smaller. From the perspective of cultural history, this dynamic environmental development played an important role in the formation of the ancient Egyptian state around 5050 cal BP (Pennigton et al., 2017).

3 Data and methods

The conceptual framework of the approach is illustrated in Fig. 2. The following sections provide information on the datasets, processing, and analysis.

Figure 2 Conceptual framework of the data processing, feature selection, and analyses. Note: NDVI signifies Normalized Difference Vegetation Index, NDWI signifies Normalized Difference Water Index, and TCB signifies brightness, TCG greenness, and TCW wetness of the tasseled cap transformation.

3.1 Data

3.1.1 TanDEM-X digital elevation model

The DEM of the TanDEM-X mission was made available by the German Aerospace Center (DLR) (see Acknowledgments) for the entire Nile Delta (coverage of approx. 260 km×180 km) (Fig. 1a). The elevation values are stored with single precision (i.e., submeter accuracy), and the vertical system refers to the height in meters above the WGS1984 ellipsoid. Note that the TanDEM-X DEM is a surface model. It therefore does not show the actual height of the terrain but the height of the land surface (including vegetation, buildings, etc.). The DEM was resampled to a resolution of 13 m×13 m.

3.1.2 Landsat archive

The Landsat Mission provides the longest remote sensing archive of optical multispectral data with the earliest acquisitions dating back to the early 1970s. In this study, we used the imagery of Landsat-5, Landsat-7, and Landsat-8 obtained between 1985 and 2019. The (passive) optical sensors acquire multispectral information on the earth's surface and offer comparable temporal, spatial, and spectral resolutions (Figs. 1b and 3).

3.1.3 Reference data

The TanDEM-X DEM and the features processed using the Landsat archive (see below) were compared to historical topographic maps of the Survey of Egypt (SoE). The maps were surveyed between 1897 and 1911 and display the topographic setting of the Nile Delta at a scale of 1:50 000. Besides the typical topographic elements (e.g., location of roads, railways, canals), information on elevated areas inside the delta is also provided by a unique cartographic signature. From the maps, it becomes clear that this information is in relation to mounts (such as, for example, tells and geziras). However, a definition of this cartographic signature is missing.

Figure 3 Study area and the number of Landsat images used in the analyses: (a) climate graph of Cairo, (b) potential and actual evapotranspiration (EVT) from the Moderate Resolution Imaging Spectroradiometer (MODIS) for selected vegetated areas in the central delta between 2011 and 2018 using the MODIS products MOD16A2-006, and (c) the number of Landsat images per year from 1985 to 2019 for July/August (summer) and January/February (winter).

Besides the SoE maps, CORONA imagery served as a reference. This high-resolution panchromatic satellite imagery was recorded in the late 1960s and early 1970s. It is available for the entire Nile Delta via the “CORONA Atlas of the Middle East” (Casana and Cothren, 2013). The datasets used in this study were acquired in the year 1968.

Finally, the high-resolution base maps provided in the geographic information system (GIS) software ArcMap served as a reference. These base maps are a compilation of recent high-resolution satellite imagery, and information is shown as true-color composites (i.e., RGB images). Most of the base maps used were recorded by the WorldView satellites, and images were acquired between 2017 and 2018 at a spatial resolution of less than 1 m.

3.2 Methods

3.2.1 Reference data georeferencing

The CORONA and the SoE datasets were georeferenced prior to the analysis using the Universal Transverse Mercator (UTM) projection at Zone 36 North and the WGS1984 ellipsoid. This was done by selecting a sufficient number (10–20) of ground control points (GCPs) for each sheet/image. The GCPs were found by comparing the datasets to the Esri base map and by identifying matching points. The root mean square error (RMSE) for both datasets was below 30 m after completing the georeferencing using an affine transformation function. The RMSEs were, therefore, sufficiently small for a comparison with the Landsat data. Finally, the georeferenced datasets were integrated into a GIS geodatabase along with the TanDEM-X DEM and the features derived from the Landsat time series.

3.2.2 Processing in the Google Earth Engine

In order to process the remotely sensed imagery, we made use of the cloud computing capabilities of the Google Earth Engine (GEE), which allows the processing and analyzing of large geospatial datasets (Gorelick et al., 2017). Processing the Landsat data and deriving the spectral-temporal metrics (STMs) relied on a GEE processing chain developed in preliminary work (Nill et al., 2019). This includes preprocessing the surface reflectance products, such as masking clouds and cloud shadows, and deriving STMs for specified features and time intervals. STMs describe the pixelwise spectral variance by reducing the temporal dimensionality into single statistical metrics such as the standard deviation. They enable the continuous coverage of large spaces.

3.2.3 Investigated features and seasons

The analyses focused on five widely used multispectral features: the Normalized Difference Vegetation Index (NDVI), the Normalized Difference Water Index (NDWI) (Eq. 1), and the brightness (TCB), greenness (TCG), and wetness (TCW) of the tasseled cap transformation (TC) (Kauth and Thomas, 1976).

$$\begin{array}{}\text{(1)}& \text{NDWI}=\left(\text{NIR}-\text{SWIR1}\right)/(\text{NIR}+\text{SWIR1})\end{array}$$

The features were processed for July/August and January/February to represent the summer and winter seasons of the year. Figure 3 shows the data availability per season and year. The number of cloud-free images available per pixel for the summer and winter seasons is shown in Fig. 1c and d, respectively.

The split into winter and summer was intended to account for different situations in terms of water supply and potential water stress due to the large seasonal difference in potential and actual evapotranspiration (EVT) (Fig. 3b). EVT was estimated from the time series of the Moderate Resolution Imaging Spectroradiometer (MODIS) for selected vegetated areas in the central delta between 2011 and 2018. The processing of the MODIS products (MOD16A2-006) was carried out via the cloud-based processing service AppEEARS (Application for Extracting and Exploring Analysis Ready Samples), which is provided by the United States Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA), among others.

Water stress is partially caused by the grain size composition of the substrate, in which infiltration rates are higher for coarser-grained deposits. Some of the features (e.g., NDWI, TCW) are known to be sensitive to the soil/plant moisture (e.g., Yan et al., 2014). As a result, locations with meaningfully different near-surface grain size compositions might be revealed through the different behavior of the STMs in winter and summer in the long run. Accordingly, the analyses and the presented results focus on the long-term differences (Δ) of the features for winter and summer, e.g., on the ΔNDWI (Eq. 2):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }\text{NDWI}={\text{NDWI}}_{\mathrm{Summer}\left(\mathrm{JA}\right)}-{\text{NDWI}}_{\mathrm{Winter}\left(\mathrm{JF}\right)}.\end{array}$$

The long-term differences were investigated for the five features listed above and for the following STMs: minimum, 5 % percentile, 10 % percentile, median, and maximum. These were processed in the temporal dimension of the stack for each feature (NDVI, NDWI, TCB, TCG, and TCW) and season (winter and summer). Therefore, the analysis investigated and compared five STMs of five multispectral features for two seasons for the entire Nile Delta using Landsat imagery acquired between 1985 and 2019.

3.2.4 Detection of anomalies

The investigations started with a proof of concept. Processing results of the GEE were compared to the results published by Ginau et al. (2019) on a buried/abandoned Nile branch in the northwestern delta close to the sites of Buto and Kom el-Arab. By using the location of the abandoned channel as a reference, all the processed Landsat features were analyzed. If and how the buried channel is revealed in the Landsat time series were then checked. This analysis aimed to select the feature that best indicates the channel in order to simplify and accelerate analysis on a delta-wide scale. Using the identified feature, anomalies (see below) were detected by visual inspection and by comparing the feature to the TanDEM-X DEM, the SoE maps, the CORONA imagery, and the Esri base maps. For all features that display the difference between summer and winter index values, e.g., the ΔNDWI, anomalies were defined by the mean and the standard deviation. In this definition, a feature value was assigned as anomalous if its value was significantly different from the mean. This significant difference, in turn, was defined by the range of 2 standard deviations centered on the mean value. For example, for ΔNDWI, the 2 standard deviation range is from −1.1 to +1.1. This means that values outside this range were defined as anomalies and were therefore of special interest. To support the visual analysis, the different images were displayed with a color bar from blue to white to red. Pixels are thus only colored (i.e., not white) if their value is outside the above-mentioned range (i.e., bluish colors for values $<-\mathrm{1.1}$ and reddish colors for values $>+\mathrm{1.1}$ for the ΔNDWI).

4 Results

4.1 Proof of concept

The processed Landsat features of both seasons are compared to the known location of the abandoned Nile branch proposed and described by Ginau et al. (2019). The branch is located several kilometers north of Buto and appears visibly in the TanDEM-X DEM as its levee is slightly higher than the floodplain (approx. 2 m) (Fig. 4). Coming from the south, the branch passes Kom el-Arab directly to the east and continues northwards towards Kom Alawi. There is no indication of this branch in the SoE maps, the CORONA imagery, or modern satellite imagery (i.e., it is just visible in the TanDEM-X data due to the aforementioned difference in elevation). At first glance, the median Landsat STMs (NDVI, NDWI, TCB, TCG, TCW) for this location do not provide information on the location of the channel. For TCB, TCG, and TCW, neither the summer nor the winter features display noticeable anomalies on the known course of the channel.

Figure 4 Comparison of the Landsat time-series features to a known abandoned/buried river branch northeast of Buto (location according to Ginau et al., 2019): (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b), (c–g) minimum (Min.), lower 5 % percentile, lower 10 % percentile, median (Med.), and maximum (Max.) of the summer NDWI values of the Landsat time series (1985–2019), (h–l) Min., lower 5 % percentile, lower 10 % percentile, Med., and Max. of the winter NDWI values, (m–q) difference (Δ) between summer and winter STM values, (r) NDWI values of the sample points P1, P2, and P3 for all Landsat images (n=825), (s) histograms of P1 and P2 for summer (n=123), and (t) histograms of P1 and P2 for winter (n=80); the ΔNDWI values are P1 $=-\mathrm{0.19}$ and P2 $=-\mathrm{0.12}$ accordingly. TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

However, the channel location is visible via lower index values in the 5 % quantile (Fig. 4d), the 10 % quantile (Fig. 4e), and, to a lesser degree, in the median (Fig. 4f) of the summer NDWI. The difference between summer (Fig. 4c–g) and winter (Fig. 4h–l) NDWI features reveals that values of ΔNDWI (Fig. 4m–q) are strongly negative (i.e., index values are lower in summer) over the channel locations, showing significant differences of $<-\mathrm{0.2}$. Ambiguities exist (especially towards the east where the SoE maps indicate a former inland water body), but the pattern of the ΔNDWI anomalies matches the proposed location and extent of the channel. Among the differences between the respective summer and winter NDWI STMs (Fig. 4m–q), the channel is best visible as an anomaly in the difference between the NDWI medians of summer and winter (Fig. 4p). This feature provides fewer ambiguities than the different images calculated using the 5 % (Fig. 4d) or the 10 % (Fig. 4e) percentile features.

The complete time series of NDWI values is exemplarily displayed in Fig. 4r for three selected sample points around Kom el-Arab. Two points lie in agricultural fields over (P1) and outside (P2) the proposed course of the abandoned Nile branch and close to each other (Fig. 4q). The third point (P3) is located in the center of Kom el-Arab. For this region, 825 Landsat acquisitions are available from between 1985 and 2019. The respective NDWI values are rather stable over time for P3 but vary along with the vegetation phenology over the agricultural fields (P1 and P2). The time series show that both fields underwent a change in land use (e.g., in irrigation or plantation) between 2012 and 2013 as the NDWI signature changes.

A comparison of the histograms for P1 and P2 for the summer (Fig. 4s) and winter (Fig. 4t) seasons shows that, in the long run, lower NDWI values (i.e., indicating lower soil/plant moisture) are more frequently found over P1 than over P2 for the summer season. For both P1 and P2, NDWI values are similar for the winter season. The difference in the seasonal NDWI median values results in ΔNDWI vales of −0.19 (P1), −0.12 (P2), and −0.05 (P3).

A similar situation is found approx. 15 km to the east close to Kom el-Khawaled (Fig. 5). For this region, another levee of a former river branch is visible in the TanDEM-X data (Fig. 5a) connecting Kom el-Khawaled with Kom el-Garad and Kom el-Nashwia (Fig. 5c).

The levee is slightly more elevated than the floodplain. Again, median summer NDWI values (Fig. 5d) are found to be lower over the levee than the winter NDWI (Fig. 5e) values. This difference is displayed as ΔNDWI in Fig. 5f. Anomalies over the levee are significant with values of −0.15.

The ΔNDWI is further compared to the general ancient (approx. 4000 BCE) delta landscape/riverscape proposed by Butzer (1976). Figure 6 shows the ΔNDWI along with a georeferenced and digitized version of Butzer's map. The mismatch in spatial resolution does not allow a comparison of both datasets in detail. However, the upper central Nile Delta (towards the apex) is indicated by strong negative anomalies of the ΔNDWI, which somewhat matches the proposed location of “sands at or near the surface” in Butzer's map, for example, the locations between Tanta and Cairo, between the modern course of the Rosetta and the Damietta branches, respectively.

The border between strongly negative ΔNDWI values and values without significant difference (ΔNDWI around zero) is, for some locations, remarkably sharp, for example, between Tanta and Zagazig, where the Damietta branch marks a clear border in the ΔNDWI values. Most of the significant ΔNDWI anomalies are negative, indicating lower summer NDWI index values. For the inner fringes of the delta and for the desert, the ΔNDWI values are close to zero (i.e., no difference in NDWI between summer and winter is observed). Positive anomalies are found less frequently but exist, e.g., for larger patterns towards Lake Burullus and Lake Manzala.

Figure 5 Comparison of the Landsat time-series features to a potential abandoned/buried river branch (location approx. 15 km east of Fig. 4) connecting Kom el Khawaled, Kom el Garad, and Kom el Nashwia (names according to the SoE Map): (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b), (c) topographic map of the Survey of Egypt (SoE), (d) median (Med.) of the summer NDWI values of the Landsat time series (1985–2019), (e) median (Med.) of the winter NDWI values, and (f) difference (ΔNDWI) between summer and winter median NDWI values. TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

Figure 6 Comparison of Landsat time-series features to the general ancient delta landscape: (a) difference (Δ) between median summer Normalized Differences Water Index (NDWI) and median winter NDWI and (b) general landscape units of the Nile Delta around 4000 BCE according to Butzer (1976). TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

4.2 Geziras and tells of the eastern Nile Delta

For the eastern delta, the region of interest lies between Bubastis (south) and Tanis (north) and between the northeastern desert margin (east) and the Damietta branch (west). Strongly negative local ΔNDWI anomalies are found in the direct vicinity of almost all larger tells and geziras in the eastern delta. This is exemplarily illustrated in detail for Geziret Sineita in Fig. 7 and for other locations in Figs. 8 and 9.

Figure 7 Processing results for the eastern Nile Delta for the example of Geziret Sineita (Tell es-Sunayta – SCA Sharqiya register 13120050, Sunaytah; EES 566 – but Geziret Bayud on SoE map – van den Brink, 1986, pp. 21, 24): (a) TanDEM-X digital elevation model (DEM), (b) Esri base map (2017/18), (c) difference (ΔNDWI) between summer and winter median NDWI values, (d) topographic map of the Survey of Egypt (SoE), and (e) values of summer and winter median NDWIs, ΔNDWI, and elevation of the DEM along profile AB. TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

Figure 8 Processing results for the eastern Nile Delta: (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b) and extent of the anomalies, (d–g) ΔNDWI (see Fig. 7) (1985–2019), (h–l) georeferenced maps of the Survey of Egypt (SoE) (ca. 1906–1912), (h–l) CORONA imagery (1968), (p–s) Esri base map (2017/18), and (t–w) TanDEM-X DEM (2012). Dashed black and yellow lines indicate the extent of ΔNDWI anomalies. Locations are, from left to right and including references, (i) Geziret Sineita, Tell es-Sunayta (SCA Sharqiya register 13120050, Sunaytah, EES 566) and Geziret Bayud on SoE map (van den Brink, 1986, pp. 21, 24), (ii) Geziret Umm Igrim/Agram (EES 569) and just south lies Tell Ibrahim Awad (SCA Sharqiya register 13120027, EES 535) (van den Brink, 1986, pp. 22, 24; van Haarlem, 2000, pp. 13–16; van Haarlem, 2019), (iii) Geziret Abu Qeih, and (iv) Tell Gumayama/Tell Gumaiyima (SCA Sharqiya register 13060032, EES 186) (Petrie, 1888, pp. 37–44; van den Brink, 1986, pp. 22, 24; Ashmawy, 2006, pp. 55–64). TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

Figure 9 Processing results for the eastern Nile Delta: (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b) and extent of the anomalies, (d–g) ΔNDWI (see Fig. 7) (1985–2019), (h–l) georeferenced maps of the Survey of Egypt (SoE) (ca. 1906–1912), (h–l) CORONA imagery (1968), (p–s) Esri base map (2017/18), and (t–w) TanDEM-X DEM (2012). Dashed black and yellow lines indicate the extent of ΔNDWI anomalies. Locations are, from left to right and including references, (i) no name, (ii) no name, (iii) Geziret Ziwilin, 1 km south of Tell Geziret Zuwelen/Tell Zuwelen/Tell Sueilin (SCA Sharqiya register 13060041, EES 324) (Petrie, 1885, pp. 29. Sect. 36; Griffith in Petrie, 1888, p. 46; Yoyotte, 1987, p. 107, n. 1; Favard-Meeks, 1999, p. 88, Fig. 1), and (iv) Gezira Samana (Habachi, 1954, pp. 479–489; Franzmeier, 2010, pp. 29–32). TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

The ΔNDWI anomalies exhibit a clear and sharp border and differences between summer and winter NDWIs are $<-\mathrm{0.2}$. Comparisons to the SoE maps, the CORONA imagery, and the modern satellite imagery indicate that these anomalies are not caused by differences in land cover (which would be visible in the modern high-resolution satellite imagery). These anomalies are also not related to former tell or gezira borders indicated in the SoE maps (i.e., almost all anomalies are found at some distance from the borders indicated in the SoE maps).

In addition, for almost all locations, the TanDEM-X DEM shows that anomalies are not caused by remarkable differences in elevation, but the terrain between the borders of the anomalies and the settlement, tell, or gezira borders is rather flat. This is exemplarily displayed in the profile lines in Fig. 7e for Geziret Sineita.

For Tell es-Sunayta (Fig. 8d), Geziret Umm Igrim, Tell Ibrahim Awad (Fig. 8e), and Geziret Abu Qeih (Fig. 8f), the ΔNDWI anomalies are asymmetric. Here anomalies are not found on both sides of the tell or gezira (i.e., one side of the tell is without a significant ΔNDWI anomaly).

The ΔNDWI anomalies are also found at smaller tells, which are indicated in the SoE maps but seem to have vanished over time. Figure 9d and e show two examples of small tells (diameters less than 250 m) in the eastern Nile Delta. These tells are indicated in the SoE maps but are not, or hardly, visible in the CORONA and the modern satellite imagery. They are, however, clearly displayed in the ΔNDWI by strong negative anomalies with values $<-\mathrm{0.15}$.

The ΔNDWI anomalies, in some cases, even encompass two or more tells that the SoE maps indicate are separated. This is, for example, the case for Geziret Ziwilin and Tell Geziret Zuwelen (Fig. 9f) and also for Geziret Umm Igrim and Tell Ibrahim Awad (Fig. 8e). For more southward locations (towards the desert margin), the patterns of the anomalies are less clear, and generally fewer anomalies of smaller extent are found. An example is the Gezira Samana (Fig. 9g). Here the ΔNDWI anomalies mostly follow the border indicated in the SoE map (Fig. 9k). There also seems to be a relation to the topographic setting at this site (Fig. 9w) as the border of the anomaly follows the terrain discontinuities.

4.3 Tells of the western Nile Delta

While ΔNDWI anomalies appear in the immediate vicinity of most of the tells and geziras of the eastern delta, this is not the case for the tells of the western delta. In fact, no significant ΔNDWI anomalies are found for the tells indicated in the SoE maps. The investigated region is located north of Damanhur, several kilometers to the east and west of the Rosetta Branch.

For example, Fig. 10 shows the ΔNDWI along with the reference data for Kom el-Ghoraf/Ghuraf (Fig. 10d), Kom el-Sheikh Ismail/Kom Abu Ismail/Kom Ismail (Fig. 10e), Kom el-Nisf (Fig. 10f), and Buto (Tell el-Fara'in) (Fig. 10g). No ΔNDWI anomalies were found that are comparable in size or shape to those of the eastern Nile Delta, although the tells are comparable in size and shape.

Figure 10 Processing results for the western Nile Delta: (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b), (d–g) ΔNDWI (see Fig. 7) (1985–2019), (h–l) georeferenced maps of the Survey of Egypt (SoE) (ca. 1897–1912), (h–l) CORONA imagery (1968), (p–s) Esri base map (2017/18), (t–w) TanDEM-X DEM (2012). Locations are, from left to right and including references, (i) Kom el-Ghoraf/Ghuraf (SCA Beheira register 100164, EES 229) (Lanna, 2005, pp. 339–363; Sist, 2006, pp. 243–249; Wilson and Grigoropoulos, 2009, pp. 173–175; Kenawi, 2014, pp. 97–99), (ii) Kom el-Sheikh Ismail/Kom Abu Ismail/Kom Ismail (SCA Kafr es-Sheikh register 090168, as Kom Ismail, EES 251) (Ballet and von der Way, 1993, pp. 12–13; Wilson and Grigoropoulos, 2009, 215–218, no. 57; Ginau et al., 2019, Fig.  5), (iii) Kom el-Nisf (EES 2), and (iv) Buto (Tell el-Fara'in, SCA Kafr es-Sheikh register 090134, EES 4), adjacent in the north to Kom el-Dahab (SCA Kafr es-Sheikh register 090171) (Hartung et al., 2009, pp. 172–188; Ballet et al., 2011; Seton-Williams, 1969, pp. 5–22). TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

4.4 Riverscape of the eastern delta

Focusing on potential former river courses, several linear and meandering anomalies are found in the eastern Nile Delta in closer proximity to Bubastis and Zagazig (Fig. 11). Local ΔNDWI anomalies are present to the west of the modern settlement of Dyarb Negm (Fig. 11c). The ΔNDWI anomalies, approx. $<-\mathrm{0.12}$ and therefore of lesser magnitude than the anomalies found near the geziras, indicate several meanders that overall stretch from south-southwest to north-northeast (Fig. 11b).

Figure 11 Processing results for the eastern Nile Delta: (a) TanDEM-X digital elevation model (DEM), (b) Landsat-8 true-color composite (RGB) (see Fig. 1b) and extent of the anomalies, (c, e, g) ΔNDWI (see Fig. 7), and (d, f, h) TanDEM-X DEM. Dashed black and yellow lines indicate the extent of ΔNDWI anomalies. TanDEM-X DEM courtesy of the German Aerospace Center (© DLR).

For this location, the orientation of the field boundaries and the terrain (Fig. 11d) also give some evidence of a former branch of the Nile. A second linear ΔNDWI anomaly is found between Zagazig and Dyarb Negm, approx. 8 km west of the modern settlement of Hihya. The anomaly stretches from the southwest to the northeast (Fig. 11b), and anomalies are strong with a magnitude of $<-\mathrm{0.15}$. None of the reference data provide further information on how this anomaly can be explained (i.e., neither the orientation of the field boundaries nor the TanDEM-X data indicate the potential location of the branch).

Finally, a third linear anomaly is found in the ΔNDWI data. It runs south of Zagazig, directly passing Bubastis, stretches from the southwest to the northeast, and therefore runs parallel to the lower reach of the more northern linear ΔNDWI anomaly (Fig. 11b). The ΔNDWI values are mostly significant, but they are of a lower magnitude with values $<-\mathrm{0.12}$. From the visual interpretation, it seems clear that there is a linear anomaly. However, the connectivity of abnormal pixels is lower compared to the other two examples, and the course of the anomaly is, therefore, more ambiguous/speculative.

5 Discussion

5.1 Detection of anomalies via Landsat time series

The detection of anomalies relied on the analysis of the ΔNDWI calculated as the difference between the long-term median summer and median winter NDWI. In the analysis the features of the tasseled cap transformation (TCB, TCG, TCW) and their differences between winter and summer seasons were also investigated. However, these features were less suited to indicate the former location of the river channels. This may result from the fact that seasonal differences between tasseled cap features primarily display differences in the surface reflection caused by different lighting conditions (i.e., reflection in winter is generally reduced as the sun elevation is lower). In contrast, NDWI is a normalized (dimensionless) index. It is therefore not sensitive to differences caused by seasonally different lighting conditions, which makes the index comparable in all the seasons.

The differencing of winter and summer features in turn aims to investigate seasonal differences. It is therefore suited to cancel out land-cover classes with temporally stable index properties, e.g., such as permanent water bodies, (dry) desert surfaces, and urban areas. This limits the analysis to land-cover classes showing seasonal dynamics due to, for example, the phenological development of vegetation or to differing moisture conditions of the plants and/or the soil. Further, the interpretation of the seasonal differences must take into account the actual land cover and its spatial variability as some land-cover types inherently cause large seasonal differences in the NDWI (see below). The approach is therefore most promising/feasible over areas displaying the same coverage.

Anomalies were identified if ΔNDWI values were outside the 2 standard deviation range centered on the mean. Strong and connected ΔNDWI anomalies were found over the known location of the abandoned Nile branch proposed by Ginau et al. (2019). It is therefore likely that Landsat features are suitable for identifying surface anomalies (over vegetated areas of the same coverage) that may be related to subsurface anomalies. In this context, the assessment indicated that ΔNDWI anomalies result, in the long run, from lower summer NDWI values. This was revealed by comparing the NDWI signatures of locations over and beside the proposed channel location. This observation is most likely caused by lower soil/plant moisture. NDWI is sensitive to these parameters (Yan et al., 2014), and summer months are characterized by the highest evapotranspiration rates, which in turn make situations with water stress or reduced water supply more likely.

Following this line of argument, and considering that anomalies are not related to field boundaries or obvious differences in land cover, they must be caused by differences in the near-surface deposits and/or in the local morphology. These are likely, as the former levee offers better drainage due to its slightly higher elevation than the floodplain. Further, the deposits are dependent on the flow velocity. They may therefore be coarser near the middle of the channel, while finer-grained deposits are expected in a more lateral position to the main channel (i.e., the surrounding floodplain) (Brown, 1997).

It is important to note that the observed seasonal differences in NDWI only become visible by investigating long timescales. A situation with reduced water supply or water stress is rather seldom for the (outer) delta also due to the irrigation. A long observation period and many acquisitions are therefore necessary to capture rather small differences in the NDWI between the two seasons. Thus, the course of the whole channel is not visible in a single image (i.e., a single NDWI), but it appears when all scenes of the time series are investigated and if summer and winter NDWIs are differenced.

However, comparisons of the ΔNDWI with the modern high-resolution satellite imagery of the Esri base map reveal that not all of the ΔNDWI anomalies can be explained by the model outlined above. There are several typical sources of ambiguities and their interpretation. Thus, ΔNDWI anomalies are frequently observed near small linear structures (e.g., roads, railroads, field boundaries) and at borders between different land-cover or land-use units (e.g., at the borders of cities). Such differences are believed to be caused by the spatial resolution of the Landsat system causing pixels with mixed coverage and by slightly different locations of the pixel footprints as the grid cells of the Landsat acquisitions are not necessarily congruent throughout the time series. Furthermore, some land-cover types and crop management practices will inherently cause large seasonal differences in the NDWI. This can be seen for plantations with broad-leaved woody vegetation, as far as this is detectable in the modern high-resolution imagery.

5.2 Anomalies of tells and geziras

Despite the limitations listed above, strong ΔNDWI anomalies were found in the surroundings of the tells and geziras of the eastern delta. The comparison to the Esri base map made clear that these are not caused by differences in the land cover as no obvious difference in coverage or usage was detected for these locations.

Considering the proposed evolution of the eastern Nile Delta (Andres and Wunderlich, 1991) it seems likely that these anomalies are related to the shallow subsurface continuation of the geziras and/or to the displaced tell and gezira material (e.g., investigated by Ginau et al., 2017, for the western delta also using Landsat imagery). Therefore, it is more likely that coarser-grained deposits are present near the surface, leading to better drainage and, on average, to lower NDWI index values in summer, as manifested by negative ΔNDWI anomalies. This is further supported as similar ΔNDWI anomalies are found over small tells that have vanished but were mapped in the SoE maps. These locations indicate that the ΔNDWI is, to some degree, sensitive to differences in the surficial substrate (Kalayci et al., 2019).

It was further observed that the ΔNDWI anomalies over the geziras in the eastern delta were asymmetric, which matches the general concept of their morphology and genesis, e.g., as demonstrated by El Beialy et al. (2001) (for Tell Tukh el-Qaramus), Andres and Wunderlich (1991) (schematic cross section for the western and eastern delta), and van Wesemael and Dirksz (1988) (for a schematic cross section of the eastern delta). This interpretation is also supported by the fact that no anomalies were found in the western delta that were similar in shape and magnitude. The tells of the western delta mostly have a different genesis and morphological setting (Andres and Wunderlich, 1991; Trampier, 2014). There are usually no large sandy bodies (i.e., geziras) or coarser-grained deposits, at least near the surface. In this context, Andres and Wunderlich (1991, p. 128) point out that fine-grained deposits in the eastern delta are embedded in a “[...] Pleistocene sand relief”, while for the western delta “[...] Pleistocene sand is not found at the surface.” While this conceptual model holds for the sites investigated by Andreas and Wunderlich (1991), later research by Wunderlich and Ginau (2016) indicated larger near-surface sand deposits at the site of Buto even though the authors note that these formations were covered by younger deposits.

Still, it is to be expected that tell material was displaced over time during the active periods (Ginau et al., 2017). Therefore, the absence/presence of ΔNDWI anomalies might finally be explained by the absence/presence of coarser-grained deposits at or near the surface whether the anomalies are due to rearranged tell material or an underground continuation of the landform.

Notwithstanding the above, the location and shape of the ΔNDWI anomalies provide interesting results for some of the tells. For example, the results suggest a linkage between Geziret Umm Igrim/Agram (EES 569) and the southern Tell Ibrahim Awad. Further, anomalies indicate a potential linkage between Geziret Ziwilin and Tell Geziret Zuwelen/Tell Zuwelen/Tell Sueilin. Favard-Meeks (1999) described the relationship of the two adjacent sites as a “double tell” with functional differences. There is a smaller settlement area in the north at Tell Geziret Zuwelen/Tell Zuwelen/Tell Sueilin and a very large cemetery area in the south at Geziret Ziwilin. The new results suggest that the two areas were connected. The northern site is much larger than previously thought, making both areas (settlement and funerary) roughly the same size. It has been observed for gezira sites that cemetery areas are located on the higher parts of geziras, while settlements are placed on lower sections (van den Brink, 1986). Possibly this is the case here, with large parts of the gezira sloping down towards the north. This may explain the reduced size in the north, which is either due to the easier leveling of this lower area or the lower parts continuing under the modern surface.

In this context, it is also interesting to note that ΔNDWI anomalies are of a smaller extent and magnitude for Gezira Samana. As indicated by Franzmeier (2010), this gezira is likely of different origin (e.g., compared to Geziret Sineita, Geziret Umm Igrim/Agram, or Geziret Abu Qeih). Surficial deposits are characterized by thin, shallower sandy deposits but more fine-grained material.

5.3 Anomalies of potential former river courses

For the western delta, the ΔNDWI anomalies were found over two abandoned river channels in the vicinity of Buto and Kom el-Khawaled. Both anomalies followed the courses of former levees. These courses were delineated using the TanDEM-X DEM as the levees are slightly more elevated than the floodplains. The observed ΔNDWI anomalies were caused by lower NDWI index values in summer, which indicate lower soil and/or plant moisture on average. This interpretation is supported, as the levees offer better drainage and are likely composed of coarser-grained deposits (Brown, 1997). Both sites display area-wide homogeneous land cover (as far as this is identifiable in modern high-resolution satellite imagery), which is an ideal setting for the detection of anomalies using the proposed approach.

The interpretation of the ΔNDWI anomalies suggests a connection between Kom el-Khawaled, Kom el-Garad, and Kom el-Nashwia and supports the proposed connection of Kom el-Arab and Kom-Alawi by a former river branch.

Similar linear, stretched ΔNDWI anomalies were found in the eastern Nile Delta in the direct vicinity of the ancient site of Bubastis and also north of Zagazig. These were interpreted as former branches of the Nile for this region, potentially of the Pelusiac or Tanitic Nile (Bietak, 1975; Tronchère et al., 2012). Compared to the anomalies found over the geziras, these were, however, less distinct (i.e., ΔNDWI anomalies were of lower magnitude), and no indication of former levees was found in the TanDEM-X DEM. The results of these anomalies therefore require more careful interpretation. They might indicate former branches of the Nile; however, the Landsat data do not allow further classification. Nevertheless, taking the information of the SoE maps into account, the potential channels have been abandoned for at least 100 years. The anomalies might indicate promising locations for future ground-based investigations, e.g., for drillings and electric resistivity surveying, which have been successfully applied to detect near-surface fluvial deposits across the delta on varying spatial scales (e.g., Andres and Wunderlich, 1991; El Gamili et al., 2001; El-Qady et al., 2011; Lange-Athinodorou et al., 2019).

6 Summary and conclusion

The paper investigated spatial-temporal metrics (STMs) of the Landsat system using remote sensing time-series data between 1985 and 2019 acquired over the Nile Delta. Comparisons to a formerly verified abandoned river branch in the northwestern delta close to Buto made clear that significant anomalies are present for some STMs over the proposed location of the former channel. The extent and location of this channel were best revealed by differencing the median NDWI values of the summer (July/August) and winter (January/February) seasons (ΔNDWI). The observed difference is likely due to lower soil/plant moisture during summer (i.e., months with high evapotranspiration). This may be caused by coarser-grained deposits and/or the morphology of the former levee. Further analyses of the ΔNDWI over the eastern Nile Delta showed similar anomalies (i.e., significantly lower summer NDWI values) in the immediate surroundings of several geziras and tells. This allowed, at least to some extent, the identification and mapping of the potential near-surface continuation of these landforms. Such anomalies were not observed for the surroundings of tells of the western Nile Delta. Additional, linear and meandering ΔNDWI anomalies were found in the eastern Nile Delta in the immediate surroundings of Bubastis, as well as several kilometers north of Zagazig. These anomalies might indicate former courses of Nile river branches (i.e., most likely of the Pelusiac and/or the Tanitic Nile). However, the ΔNDWI does not allow an unambiguous interpretation. The analyses have shown that the Landsat archive is promising for (geo-)archeological questions, especially in the context of landscape archeology. A transfer of the methodology to similar environments seems feasible considering the relatively simple approach, the global availability of the Landsat data, and the benefits that arise from the cloud processing in the Google Earth Engine. Future research will continue analyzing the ΔNDWI for the entire delta but will also investigate other spatiotemporal metrics that can be deduced from the Landsat time series (cf. Ginau et al., 2017). Furthermore, the benefit of including additional imagery from passive (e.g., MODIS, Sentinel-2) or active (e.g., Sentinel-1) remote sensing sensors will be assessed.