Use of Time-Series NDWI to Monitor Emerging Archaeological Sites: Case Studies from Iraqi Artiﬁcial Reservoirs

: Over the last 50 years, countries across North Africa and the Middle East have seen a signiﬁcant increase in dam construction which, notwithstanding their beneﬁts, have endangered archaeological heritage. Archaeological surveys and salvage excavations have been carried out in threatened areas in the past, but the formation of reservoirs often resulted in the permanent loss of archaeological data. However, in 2018, a sharp fall in the water level of the Mosul Dam reservoir led to the emersion of the archaeological site of Kemune and allowed for its brief and targeted investigation. Reservoir water level change is not unique to the Mosul Dam, but it is a phenomenon affecting most of the artiﬁcial lakes of present-day Iraq. However, to know in advance which sites will be exposed due to a decrease in water level can be a challenging task, especially without any previous knowledge, ﬁeld investigation, or high-resolution satellite image. Nonetheless, by using time-series medium-resolution satellite images, combined to obtain spectral indexes for different years, it is possible to monitor “patterns” of emerging archaeological sites from three major Iraqi reservoirs: Mosul, Haditha and Hamrin lake. The Normalised Difference Water Index (NDWI), generated from annual composites of Landsat and Sentinel-2 images, allow us to distinguish between water bodies and other land surfaces. When coupled with a pixel analysis of each image, the index can provide a mean for highlighting whether an archaeological site is submerged or not. Moreover, using a zonal histogram algorithm in QGIS over polygon shapeﬁles that represent a site surface, it is possible to assess the area of a site that has been exposed over time. The same analyses were carried out on monthly composites for the year 2018, to assess the impact of monthly variation of the water level on the archaeological sites. The results from both analyses have been visually evaluated using medium-resolution true colour images for speciﬁc years and locations and with 3 m resolution Planetscope images for 2018. Understanding emersion “patterns” of known archaeological sites provides a useful tool for targeted rescue excavation, while also expanding the knowledge of the post-ﬂooding impact on cultural heritage in the regions under study.


Introduction
Among the many threats to which archaeological sites are subjected in the Middle East and North Africa (MENA) region, dam constructions and artificial lake formations have been some of the most common over the last 50 years (Figure 1). During this time, hundreds of dams have been built in the Middle East, providing undeniable benefits but also having an impact at the socio-economic, environmental, and geopolitical level [1][2][3], [4] (pp. [13][14], [5] (pp. [4][5][6]. Archaeological sites were also impacted by the construction of dams and subsequent formation of artificial lakes, a threat soon recognised by national and international institutions. Several salvage projects, consisting of archaeological surveys and rescue excavations, have been carried out in different regions where cultural heritage was endangered by the construction of dams [6,7].  [8,9]. Some of the best known and best-published examples for Near Eastern archaeology are the salvage surveys and excavations carried out from the late 1960s until 2010-20112010s along the Euphrates and Tigris rivers (e.g., for the Euphrates, see refs [6,7,[10][11][12][13]. For the Tigris river see, e.g., refs [14][15][16][17]) ( Figure 2). However, two studies have pointed out that the general awareness of the potential loss of cultural heritage linked to dam construction has been less than optimal, even if archaeological survey projects were carried out in the past [18][19][20]. Using Turkey as a case study, Marchetti et al. [19] describe how heritage protection still lacks a systematic approach to the safeguarding of archaeological sites and features, both pre-and postflooding [19] (pp. [20][21], [21] (p. 13), [20] (pp. [30][31][32]. Regarding the Euphrates, they point out that only 12 out of the 23 areas that would be affected by dams along the river have After this brief overview, it is clear that there is evidence of a significant loss of archaeological data due to water submersion. However, we have evidence that some sites might be still accessible under specific conditions of receding waters. As explained below, these conditions can be prompted by natural or anthropic events, such as climate change or water management policies. Identifying and monitoring water bodies has been a common practice in disciplines such as geography or geology, and one of these studies can be used as a starting point for this research. A study group from the University of North Carolina explored the impact of climate-and politics-related events on some artificial lakes of present-day Iraq, with a specific focus on the Haditha and Mosul reservoirs [27]. Using remote sensing data, they explored lake-surface changes from their formation in 1985 and 1986 to 2016. The result was a better understanding of how a combination of anthropic and natural events brought drastic changes in the lake-surface areas over the last 30 years. Sudden changes, registered for the years 1991, 2001, 2009, 2015, were mostly linked to water management policies tied to political events (e.g., the series of conflicts from the Iraq-Iran war to the ISIS invasion, [27], pp. 271-273). On the other hand, prolonged years with low water level were linked to a reduction in precipitation along the Euphrates and Tigris rivers, coupled with the number of dams upstream of each reservoir. A more widespread drought was registered in the area since 1999, with a peak in 2009, which in turn affected the water level of the reservoirs [4] (pp. 58,59), [27] (pp. 273-276), [28] (pp. 270-276).
Usually, remote sensing research focused on studying submerged archaeological features relies on methods which are difficult to apply in countries such as Syria and Iraq, e.g., airborne laser remote sensing [29] or sonar technologies mounted on boats [30,31]. In the Near East, research has focused mainly on mapping relict water features [32][33][34][35][36] or used coarser-resolution imagery that does not allow for immediate validation of the results [37,38]. Instead, cultural heritage monitoring using remote sensing techniques focused mainly on looting and destruction caused by conflicts [39][40][41][42][43][44][45][46] (the EAMENA project also released an online database with the location and types of damage affecting archaeological heritage, available at http://eamena.arch.ox.ac.uk/ (last accessed on 14 December 2020). On the other hand, other studies that focused on water body detection or temporal change mapping [47][48][49][50][51][52] were never linked to archaeological heritage research, even if we now know that sites can be affected by these phenomena as well (on this topic, using the example of the Delmej reservoir in southern Iraq see Marchetti et al., 2018 [21], pp. 12-15; on temporal change detection for mapping damage to archaeological sites caused by human activities but not focused on dam construction, see also Rayne et al., 2020 [53]). For these reasons, the aims of this paper were:

•
To evaluate whether more archaeological sites than the two mentioned before are affected by the cycles of water retraction; • To establish a systematic and easy-to-reproduce methodology for the mapping and monitoring of these events over a long time span, on an annual basis. The selected time span will cover from the completion of the dams to the present day, but the method could be adapted to different situations; • To evaluate the impact of the interannual variability of the water level and its impact on the archaeological sites, by monitoring, monthly, over a single year.
The monthly analysis was carried out over the year 2018, which was chosen because a sharp lowering of the water level was registered in all the three reservoirs. The examination of the water level fluctuation and its impact on archaeological heritage is not only interesting from the research point of view, but it can also be useful in helping future targeted excavations or preservation activities.

Study Areas
Given the evidence quoted above, the choice of the Haditha and Mosul Dam lakes as study areas seemed obvious. Moreover, both areas were surveyed before the construction of dams, and archaeological reports are available for each of them (see details in Section 2.2). However, there is a third reservoir that can be included as well: the Hamrin Dam area was surveyed before its flooding in 1981, and many archaeological sites were recovered beforehand [54,55]. As described below (Section 2.1.3), the dynamics of the resulting artificial lake make the Hamrin reservoir another ideal candidate for this study.  The Haditha lake is prone to strong surface reductions: for example, during  the droughts recorded in 2007-2009 and 2013-2014, the surface lake area dropped by 72% and 62%, respectively [27] (pp. 271, 272).

Mosul Dam
The Mosul Dam area, located along the Tigris river, extends roughly from the northwest of the reservoir (37 • 0 27.9792" N, 42 • 21 5.1624" E) to the south-east of the dam (36 • 33 32.49" N, 43 • 1 57.3096" E) ( Figure 5). The dam itself, the largest in Iraq, was completed in 1986, creating the Mosul reservoir, (previously known as Lake Saddam) which covers an area of approximately 372 km 2 with a capacity of 11.2 km 3 . The Mosul Dam lake, if compared to the Haditha Dam, is less prone to substantial reductions in water, possibly because of the presence of fewer dams upstream. Periods of substantial water loss are more or less the same as the former dam, but with a significant loss in 2011, when the surface lake area dropped by 50%, caused by a water management policy aimed at avoiding a dam failure [27] (p. 276).

Hamrin Dam
The Hamrin Dam area, contrary to the other two areas, does not extend along the two main Mesopotamian rivers, but instead along the Diyala river, one of the tributaries of the Tigris. The area extends from the north-west of the reservoir (34 • 22 57.126" N, 44 • 49 16.4928" E) to its south-east (34 • 2 1.9392" N, 45 • 12 30.6576" E) ( Figure 6). The dam was completed in 1981, forming the Hamrin reservoir, which extends for approximately 340 km 2 with a capacity of 4.2 km 3 . The Hamrin lake, unlike the former two reservoirs, is subject to frequent surface area changes, sometimes affecting the very shape of the lake. Periods of low water level are similar to those of the other two dams, especially regarding the effects of widespread droughts (1999-2001; 2009; 2013-2014). However, the most substantial reduction in the lake area was registered in 2008, when it dropped by 80%, due to the damming of the Lund river in Iran [4] (pp. 139-140), [52] (pp. 58, 59).

Archaeological Sites
This research also relies on the knowledge of site locations; however, accurate coordinates for archaeological sites are not always available. This is especially true for settlements investigated by older salvage projects which, due to time constraints, rarely revisit a site after the first discovery, especially if it was not a target for rescue excavation. Moreover, technological limitations impact the quality of information about a site location [56] (pp. [47][48][49][50]. In order to work properly in a GIS environment, data on archaeological sites from the study regions were gathered from two primary sources:

•
An open, and freely available dataset of Google Earth placemarks, available as a .kmz file named "ANE.kmz", which stores name and location of more than 2500 archaeological sites across Egypt and the Near East The addition of the second source of data was necessary because the ANE.kmz, although massive, does not store a complete inventory of all the archaeological sites in the study areas.
In order to digitise sites from the distribution maps with enough accuracy, it was necessary to operate in a series of steps, which are part of a common procedure for repositioning sites in GIS, especially for Near Eastern archaeologists, who often deal with old survey data. These steps, apart from georeferencing the maps, combine the analysis of site location descriptions and topographic maps provided by scientific reports and the use of satellite imagery, in order to improve the accuracy of repositioning a site in GIS (see ref [56] pp. 50-53 for an overview of the process). However, it was not always possible to relocate all known sites in each area with enough confidence, and the resulting dataset used in this research is necessarily incomplete.
The final dataset consisted of 54 sites for the Haditha Dam, 48 for the Mosul Dam, and 42 for the Hamrin Dam, for a total of 144 archaeological sites (see also Section 2.4.2). The dataset used here for the Mosul Dam is different from the one published in Sconzo; Simi, 2020 [60]. While location precision for sites analysed in this paper was checked against that dataset, the overall accuracy and related results might have changed. The Mosul Dam dataset published in Sconzo; Simi, 2020 [60] will be the target of a future, more accurate, endeavour conducted by P. Sconzo, F. Simi, and A. Titolo using the same methodology applied here.
This research also relies on the possibility of measuring the resurfaced extent of an archaeological site. For this reason, polygons representing the extent of a site were drawn in GIS, based on their visibility on satellite images coupled with topographic maps, when available. It was not possible to reconstruct the actual extent of all the sites available; therefore, polygons were drawn only for those sites that satisfied the following conditions:

•
On satellite images, the site shows evident features such as walls that delimit its extension; • Visible topographic features (e.g., a Tell) that was taken as minimum site extent.
Such conditions were necessary in order to account for an ongoing discussion about the reliability of overall sites dimensions reconstruction and definition of clear limits for archaeological settlements. This discussion is also tied to the matter of how Near Eastern archaeological surveys have historically been conducted (see, e.g., Banning, 2002 [61], pp. 81-84; Lawrence, 2012 [56], pp. 54-56). All these elements dictate careful consideration when dealing with sites' extent. For this reason, simple dimensions such as, e.g., "50 × 100 m" provided by survey reports were not considered here, and only sites that met the conditions highlighted before were included in the analysis. It was possible to reconstruct the surface extent of 51 sites, of which 11 were for the Haditha Dam, 9 were for the Mosul Dam, and 31 were for the Hamrin Dam.

Satellite Images
The goal of this research was to monitor events over 40 years, but also to be as easy to reproduce as possible; therefore, satellite images were chosen accordingly. To cover a time span from 1984 to the present day, it was an obvious choice to rely on Landsat satellites, and specifically Landsat 5, 7 and 8. Some Landsat 5 images were also acquired as replacements for some Landsat 7 images that were too problematic to use due to the Scan Line Corrector (SLC) failure. For more recent years (from 2015 to 2019) Sentinel-2 images were acquired. The result was an average of 36 images for each area, but because it was necessary to produce Normalised Difference Water Index (NDWI) values for the year before the completion date of the Hamrin Dam (1981), one Landsat 2 image was acquired as well, only for that area ( Table 1). All images were taken at the Top of Atmosphere (TOA) reflectance (for details, see Section 2.4).
The monthly interannual analysis for 2018 was carried out on 12 Sentinel-2 images for each area, for a total of 36 images.
To obtain a visual confirmation of the results from the pixel analysis of the NDWI images, it seemed more convenient to use higher resolution imagery, at least for this last stage. A good compromise between cost-effective and good resolution imagery is provided by Planetscope images, distributed by Planet at no cost after a European Space Agency (ESA) Third Party Mission (TPM) application. The Planetscope satellite constellation is formed of more than 120 individual satellites (named DOVE) that cover the globe daily at a 3-3.7 m Ground Sampling Distance (GSD). More information is available at: https://www.planet.com/products/ (last accessed on 14 December 2020). Access to Planetscope images was granted for the entire year of 2018 and limited within the three study areas indicated above. Between specific close-up scenes and monthly snapshots, 46 Planescope Ortho Scene (Level 3B) images were gathered through the Planet Explorer portal ( Figure S1, the portal is accessible at https://www.planet.com/explorer (last accessed on 14 December 2020)). For previous years, Sentinel-2 or pansharpened Landsat 8 images were used instead.  Figure 7 illustrates the methodological approach adopted in this study. Firstly, an assessment of the processing methods useful to discriminate between submerged and emerged areas was made, and the NDWI was chosen as the most appropriate. Secondly, NDWI images were generated in Google Earth Engine, and then imported in QGIS for the pixel analysis. This latter stage was carried out using Saga's Add Raster Values to Points algorithm on the archaeological sites point shapefile. The analysis of the extent of the emerging sites was carried out in QGIS as well, and it required a reclassification of the NDWI images beforehand, paired with the use of the Zonal Histogram algorithm on the polygon shapefiles. The results obtained were then checked with Planet satellite images for the year 2018, and with Sentinel-2 or Landsat images for any previous year not covered by the assigned Planet quota. To provide reproducibility, two QGIS models and one R script (.rsx) to use in QGIS were created; moreover, the analysis was also reproduced in R, and both the R code and the QGIS models and code are linked in the Supplementary Materials.

NDWI
Medium-resolution images alone can rarely be used to visualise archaeological sites; hence, to achieve the aims of this research, there was a need for a method able to clearly discriminate between land and water surfaces.
There are different techniques available to delineate and map water features (for an overview, see refs 48,63), but one of the most widely used in is the Normalised Difference Water Index (NDWI). As argued by Du et al. ( [48], pp. 672, 673), NDWI is a suitable tool for mapping water features, because it is easier to obtain and generally more reliable than single-band classification methods, such as density slicing (on this topic, see also Qiao [64]; the same approach was also used by Marchetti et al., 2019 [19]). Both qualities perfectly fit the scope of this research, and this index was chosen as the one to delineate water masses.
Two formulae are generally used to obtain the NDWI, by combining the green band with either the near-infrared band (NIR) or the short-wave infrared (SWIR) band of satellite images. These formulas are usually called after their creators, or named NDWI [65] and modified NDWI [66], respectively: In the NDWI images, values close to or below an arbitrary threshold (usually 0, but for other examples see Ji et al., 2009 [67]) can be identified as land formations, while values above the threshold can be classified as water. There is no unanimous agreement about which of the two formulas works best for any occasion, and literature about their use mostly shows how the benefit of the two bands might be very situational and tied to the immediate surrounding environment of the water bodies [50,[68][69][70][71][72]. The formula used here is Equation (2), because it is generally agreed that, unlike Equation (1), it can better distinguish between water and built-up areas without returning false positives [66] (pp. 3026, 3027). One weakness of this formula is, however, that it tends to return positive values when snow is present [51]. The selected study areas do not have dense vegetation or snow around them but can present sparse built-up features close to the lakes; therefore, Equation (2) seemed more appropriate. The NDWI is generally aimed at identifying water bodies, but if the location of an archaeological site is known, then the index can be used to check if that location is under water or not. In fact, on satellite images, the basic difference between an emerged and a submerged archaeological site will be whether pixels can be identified as representing water or not.

Google Earth Engine
The processing to obtain the NDWI images was carried out in Google Earth Engine (hereafter GEE), because it is capable of processing a large quantity of data in a short time [73]. The number of raw images necessary would have also required significant amounts of space on any physical drive, a problem that can be avoided by using the cloud-based GEE platform. A script was written to generate the NDWI images for every year and areas using Equation (2) seen at Section 2.4.1. Each image was acquired at the TOA reflectance, and a cloud masking function was applied to the data as well. During the process, a filter for a low cloud pixel percentage was applied; however, this percentage was tweaked in cases when there were no available images that satisfied those criteria. Nonetheless, there were still some problematic years where images without dense cloud cover were available, e.g., 1985 for the Haditha Dam, and December 2018 for the Mosul Dam. After the image selection and the generation of the NDWI, the annual composite was created by taking the median of the NDWI values at each pixel through the collection of NDWI images in each year interval [71] (p. 158). The same procedure was applied for the generation of the 2018 monthly composites. One NDWI composite was also generated for the year before the construction of each dam. The Hamrin Dam was completed before the first Landsat 5 satellite acquisition (1984); therefore, the NDWI image was generated using a Landsat 2 image. Given that the SWIR band was not present on the Multispectral Scanner System (MSS) sensor of Landsat 2, the NIR band (Band 7) was used instead, and the image underwent a transformation process from digital numbers (DN) to reflectance values (the process for the Landsat 2 images was compiled on a separate script). The low cloud filter and the cloud masking function were applied to minimize potential impacts from the atmosphere [74] (see the discussion Section 4), while the median usually helps in removing noise and minimizing the effects of clouds and shadows [53]. In R, the process used in GEE was carried out using the RGEE package [75]. After obtaining all the processed images needed, the subsequent analyses were carried out in QGIS.

Pixel Analysis and Zonal Histogram
In QGIS (versions 3.10/3.16), the pixel analysis was carried out using Saga's Add Raster Values to Points algorithm. This algorithm allows for the retrieval of pixel values at each point location (i.e., the archaeological sites) across multiple selected raster images. This process was applied to the archaeological site shapefiles using the annual NDWI composites, and then repeated with the 2018 monthly composites.
The algorithm generates a new shapefile with the pixel values at the point location in all the underlying images. In the resulting shapefile, using 0 as a threshold between water and land, if a point has a pixel value below 0 in all the images, it can be assumed that the area was likely not affected by the dynamics of water retraction, while if the values are always above 0, the area was always submerged. Furthermore, if the point shows alternating values above and below 0, it can be assumed that part of the site re-emerged more than once due to the lowering of the water level. Although it can be argued that a single point is not very representative of the entire area of an archaeological site, its use can be suitable to give a first idea of the dynamics affecting an archaeological site, because it gives an unambiguous and quantifiable value. The point was usually placed at the highest point of a site, which is more prone to emerge after a change in the water level. Moreover, once a value below 0 is registered, a visual inspection of the area around the sites using satellite images can help validating the results and to assess the extent of the emersion.
The next step was to assess to what extent the sites were impacted by the water retraction, both over the long term and during the sample year of 2018. The NDWI images were reclassified so that all the pixel values between −1 and 0 were counted as 0, and all the values between 0 and 1 were counted as 1. The reclassification allowed for the easy application of the Zonal Histogram algorithm to count the number of 1 and 0 pixels within each polygon representing site extents. In R, the pixel analysis was carried out using the raster package [76], while the qgisprocess package was used for the zonal histogram analysis. This last package is still in an experimental phase, but it is available at https://github.com/paleolimbot/qgisprocess (last accessed on 25 January 2021).

Results
Results show that a larger number of sites were affected by the cycles of water retraction than just the two presented at the beginning of this paper (Figure 8). The overall number of sites for which a negative pixel value was registered at least once during the time span analysed was 73 ( Figure 9). Most of these sites were located at the edges of the three reservoirs; however, some differences can be highlighted for each area, especially regarding the number of sites that emerged per year.

Haditha Dam
Of the 54 sites analysed for this region, 27 sites registered a negative pixel value at least once from their submersion in 1985, which means that it is likely that they have been affected by the dynamics of water retraction. Most of these sites are located in the western part of the reservoir, but others were also evident further to the east ( Figure 10 and Figure S2). Following what was already known about the years with lower water level (Section 2.1.1), 2009, 2010 and 2015 had the highest number of resurfaced sites, with 25, 13 and 24 sites, respectively (Figure 11a-c). Another significant year was 2018, when 17 sites resurfaced due to another water lowering event (Figure 11d). On the contrary, other years in which the lake was known to have suffered a retraction (see Section 2.1.1) showed only a limited impact on archaeological sites: in fact, in 1991, four resurfaced sites were registered, while five resurfaced in 2001.
In general, while years with exceptionally low water levels have a high number of resurfaced sites, some archaeological sites are more prone to emerge even if the water level is not low (Figure 12). Among these, there is Kifrin, which emerged from the water 11 different times. Analysis of the resurfaced extent of the site showed that even when the water level was not extremely low, (e.g., 1999, 2001) more than 80% of the site was out of the water ( Figure 13). Other interesting results come from the site of Sur Telbis, an important Iron Age and probably Parthian site [10] (pp. 356-358), which was out of the water four times, in 2009-2010, 2015 and 2018. A high-resolution image in Bing Maps likely from 2011 shows that it suffered significant erosion on the river side, and only the highest part of it survived, probably because it was never directly reached by water ( Figure 14). This is confirmed by the analysis on the polygon shapefiles, which shows how the emerged extent of the site changed every year, suggesting once again substantial erosion damage. Another interesting result is visible for the Iron Age site of Glai'a [10] (pp. 166-180). It is located in the middle of the reservoir and was recorded as being out of water three times, in 2009, 2015 and 2018 ( Figure 15). The images from 2018 show that the south-eastern portion of the site was visible, together with its double wall. Given its position in the middle of the reservoir, the site rarely experienced significant emersions; however, in 2009, 72% of the site was out of the water (see also Section 3.4).

Mosul Dam
Compared to the Haditha Dam, in the Mosul Dam region fewer sites were affected by the water-level reduction cycles. In fact, in the Mosul Dam area, 16 sites emerged during the chosen time span, while 23 do not show any sign of resurfacing ( Figure 9). Due to the shape of the lake, the sites that emerged are more evenly distributed than in the Haditha Dam reservoir; they were found on both the eastern and western parts of the artificial lake ( Figure 16 and Figure S3). A significant number of settlements resurfaced in similar years to the Haditha Dam, for example 2009, 2011 and 2015 (with 9, 16 and 12 sites, respectively), but, as before, it is possible to add the year 2018 as well, when 16 sites emerged from the ebbing waters (Figure 17a-d).  One significant difference in respect to the Haditha reservoir is that archaeological sites at the edges of the Mosul Dam registered a negative pixel value more often than the sites in the former dam. The result is that some of them were probably above the water level more often than the time they were below it ( Figure S3). Among these sites, there are Ger Matbak, Khirbet Hatara, and Tell Baqaq 1 (26, 20, and 31 times, respectively).
For Ger Matbak, the polygon analysis showed that it was rarely completely submerged, and that even years with moderate water levels (e.g., 1998 and 2002) saw 60% of the site emerge (Figures 18 and 19). Moreover, this latter analysis showed that the entirety of the site was out of the water for four consecutive years (2014-2018). The area of Kemune site, mentioned above (Section 1), was recorded to be above the water level four times during the 35 years of analysis. It is possible that the exceptional preservation state of the site might be due to the fewer number of times it was subject to erosion by the ebbing waters (see Section 4). A note on this site is necessary, however: no coordinates were available at the time of the analysis; therefore, its position was assessed from the published photographs, available satellite images, and from the dataset published in Sconzo; Simi, 2020 [60]. It is likely that the pixel values and the respective emersion times will differ with more precise coordinates. The same can be said about other archaeological sites and the general results for the region, which might change with a more accurate and extended dataset. Other examples of emerging sites come from some settlements such as Khirbet Kharasan (9 times) and Tell Dhuwaij (12 times). In contrast to Ger Matbak, Khirbet Kharasan was most affected by significant reductions in the water levels; its resurfaced extent was around 80% only on two occasions, in 2011 and 2018 (Figures 18 and 20).

Hamrin Dam
This area shows quite different results in respect to the other two reservoirs; while the Haditha and Mosul reservoir maintained their overall shape even in years with very low water levels, the Hamrin lake was subject to substantial shape changes. For this reason, the number of sites resurfacing from the lake was very high when compared to the overall number of settlements in the region used in this analysis. Of 42 sites analysed, 12 had never been affected by artificial lake formation, while 30 of them had resurfaced at least once, meaning that all sites have resurfaced at least once ( Figure 9). During the drastic change in water level in 2008 (see Section 2.1.3) even the sites in the middle of the reservoir were above the water level. Other periods with a significant number of sites emerging from the reservoir were registered in 2000 (23 sites) and in 2015 (17 sites); in 2018, the water level still impacted the area, but was not as low as the other years, and in turn resulted in fewer sites having resurfaced (12) (Figures 21 and 22). Many sites show frequent emersion rates ( Figure S4), such as Tell Yelkhi (32 times), or the group of settlements around Tell Razuk, which seem to have been out of the water more often than they were submerged (more than 28 times). Other settlements that have resurfaced less often are, e.g., Tell Zubeidi (13 times) or the group of sites around Tell Baradan: here, Tell Hadad and Tulul es-Sib emerged 13 and 15 times, respectively, while Tell Baradan itself resurfaced 31 times. The polygon analysis showed that both Tell Yelkhi and Tell Baradan were affected by most of the fluctuations of the water level, and that they could have most of their surface exposed over one year while being almost completely submerged the next, or vice versa. This happened, for example, from 1998 to 1999, or again in 2006 and 2007; moreover, it also means that they probably suffered significant erosion damage (Figures 23-25). Regarding the other sites mentioned above, Tell Razuk rarely showed a sharp drop in the emerged extent, and Tell Hadad showed a pattern similar to that of Sur Telbis in the Haditha Dam, with a resurfaced area that changed almost every year, suggesting once again a significant erosion effect ( Figure 23).

Interannual Variability-2018 Monthly Analysis
The analysis of the 2018 monthly NDWI images showed that interannual variability can have an impact on the extent of the resurfaced sites. In the Haditha and Mosul Dam reservoirs, 2018 was a year of particularly low water levels, and both reservoirs showed a comparable behaviour in the timing of their retractions (Figure 26a,b). During this year, months of high-water level were associated with the beginning of Summer, during May and June, while the minimum levels were registered in October and November for that year. The Hamrin Dam lake followed a similar trend, even if the period of lower water level began slightly earlier, during September (Figure 26c). Individual pixel values showed a large number of sites likely out of the water (i.e., with negative pixel values) for all the three reservoirs. During October and November, the Haditha Dam registered 26 and 25 emerged sites, respectively, also counting 13 which showed a negative pixel value for the entire year. In the Mosul Dam lake, 19 sites of the 48 analysed were out of the water for both the same months, including four sites with a negative pixel value for the entire year. Unfortunately, December was a particularly cloudy month in the Mosul Dam region, which prevented the extraction of any pixel values for this month. In the Hamrin Dam region, the earlier water reduction resulted in September and October as the months that registered a great number of sites resurfaced, 17 and 20, respectively, including eight sites that could be considered as out of the water for the entire year.
In terms of spatial distribution, the visible pattern is similar to that of the long-term analysis, with sites at the edges of the reservoirs being more prone to emerge. However, given the particularly low water level in the Haditha and Mosul Dam during this year, even sites closer to the central part of the reservoirs emerged at least once (Figures 27 and 28). A similar pattern is evident for the Hamrin Dam, which showed a spatial distribution of emerged sites close to that of some previous years of low water level, such as 2015 and 2010 ( Figure 29).   The results from the point pixel analysis can be integrated with the analysis of the resurfaced extent: a pattern of decreasing and increasing emerged area was indeed registered for many sites in all three dams from May to December 2018 (Tables 2-4). In the Haditha Dam reservoir, while many sites were almost completely out of the water during the entire year, it is possible to appreciate the progressive emerging and submerging of a site such as Glai'a. This site was mostly submerged during summer, but approximately 80% of it resurfaced during October-November (Figure 30j,k). A comparable behaviour has been registered for Khirbet Kharasan in the Mosul Dam lake (>80% of site emerged) ( Figure 31). In the Hamrin Dam sites such as Tell Yelkhi, Tell Khesaran and Tell Kharbud, all followed a similar trend, with more than 85% of their extent emerged during September-November ( Figure 32). Planetscope images, acquired for specific sites for each month in 2018, confirmed the results of the analysis; the progressive emergence of the sites was clearly visible on the high-resolution images.

Accuracy Assessment
Accuracy assessment for time-series analysis is notoriously difficult [77] (see also Congalton; Green, 2019 [78], p. 233), and in this case also complicated by the lack of field recording, the less-than-optimal coverage, and the still-prohibitive costs and political restrictions for obtaining high-resolution images from commercial satellites. For these reasons, accuracy assessment was carried out on pansharpened true-colour composites from the same satellites (a similar approach due to the same limitations was adopted by Rayne et al., 2020 [53]). While this is not always optimal, because it mostly assesses the accuracy of the reclassification for that single area [79], it was the best option especially for less recent years [78] (p.235) (see also Olofsson et al., 2014 [77], p. 14) and it also provides consistency between reclassified and reference data. Pansharpening is not possible on Landsat 5; therefore, validation has been carried out on 30 m true-colour composites (thus, validation results should be taken with care). For more recent years, when possible, high-resolution satellite images from Bing Maps or Google Earth were used. For the monthly composites, pansharpened composites and Planetscope images were used as reference data [53].
Given the number of images and the different sensors used in this research, validation was carried out for two images per sensor, one with a lower water level, and one with a higher water level. Sampling was carried out at the per-pixel level using a stratified random sampling with two classes: water, and non-water pixels [77,78]. This method allows for allocating enough samples for each class depending on their area proportion, thus accounting for the lake area variation between images. For the annual images, validation was carried out on the entire research area. When mixed pixels were detected, the reference class was assigned to either water or non-water depending on the most prevalent class inside the sample. For monthly composites, validation was limited in areas covered by the Planetscope data, around 20 km maximum around some archaeological sites (Glai'a in the Haditha Dam, Khirbet Kharasan in the Mosul Dam, and Tell Yelkhi in the Hamrin Dam, see . All the validations were carried out in QGIS, using the Semi-Automatic Classification plugin [80]. Overall, the NDWI showed promising results for all the areas and satellites used (Table 5), as usually expected by a simple water/non-water classification in areas with a clear separation between a lake and its surroundings [81]. Errors of omission were generally rare and were mostly tied to the medium resolution of the images used for the NDWI or to the surroundings of each lake. In fact, agricultural fields around the Hamrin Dam lake increased in number since the last decade, and more recent Landsat or Sentinel-2 images did not identify smaller canals. However, in the context of highlighting emerged areas from the lakes, these errors had a minimum impact on the results. This is especially evident when assessing the accuracy at a smaller scale for the monthly images, in which the omission errors of the Hamrin Dam were significantly lowered. Errors of commission were slightly more frequent, but again, mostly tied to the medium resolution of the images used, sometimes unable to distinguish mixed pixels in areas of finer sediments at the edges of the lakes. It should be mentioned that no builtup pixels were incorrectly registered as water, thus confirming the correct choice of the NDWI xu for these regions (see Section 2.4.1). Nonetheless, as shown in Fisher; Danaher, 2013 [81], commission errors can be potentially lowered by experimenting with different thresholds. However, this should be evaluated for each area in order to find the respective optimum threshold, which is something beyond the scopes of this paper.

Discussion
By including more than one case study, it was possible to see how the different dynamics of each lake affect archaeological sites. The Mosul and Haditha lakes are mostly similar in these terms, while the Hamrin lake is a peculiar case, much more unstable, and with its sites that are likely to emerge more often than in the former two areas.
The results showed that the cyclical water retraction is indeed affecting a significant number of archaeological sites. This number might even be higher if we consider that older archaeological surveys suffered from visibility issues and methodological limits, and that this research had access only to a limited dataset (see Section 2.2).
The analysis shows how the sites in each area can be divided into: • Sites that emerged cyclically from the waters; • Sites that were never affected by the reservoir; • Sites that were always submerged during the observed time span.
Among the first category of sites, the one of interest for this research, the pixel analysis was able to identify those that were more likely to emerge more often. Usually, these sites are found at the edges of the reservoirs and are therefore impacted by any water level change; this is the case for sites such as Kifrin in the Haditha Dam, Khirbet Karhasan in the Mosul Dam, and Tell Zubeidi in the Hamrin Dam. On the other hand, drastic and substantial water level changes can directly affect inner sites, which can emerge, partially or entirely, from the waters even for a prolonged time. This is the case for sites such as Glai'a in the Haditha Dam, Jamrash and Kemune in the Mosul Dam, and Tell az Zawyeh in the Hamrin Dam.
There are, however, some considerations to make about the results presented above, tied to the limits of a remote sensing approach.
Firstly, even if an archaeological site is shown to be above the water level, the access to it might still be problematic, e.g., the environment might still present marshes, unstable sediment loads, or it might still be submerged. This is the case for sites such as Tell Yelkhi ( Figure 32).
The current state of an archaeological site after its emersion will also be tied to its characteristics before the formation of the lakes. For example, a site characterised by a surface scatter of archaeological material will have fewer chances to survive than a site with standing architectural elements (even if artefacts might be lost in this case as well).
The dynamics of the receding waters also have an impact on the preservation of a site. In fact, frequent emersions and submersions might damage a site more than its long submersion, through processes of erosion, particularly dangerous for sites with subtle archaeological features [22] (p. 72). This effectively means that, even if a site emerges from the water, its archaeological features might not have survived.
A visual inspection of the results is always advisable, because the pixel analysis might return negative values (i.e., land masses) for areas with high sediment loads. This is not a classification error, because the NDWI correctly identifies sediments as non-water pixels, but it is something to consider when planning a field inspection. This is evident for the Island of 'Ana, a multiperiod archaeological site in the Haditha Dam reservoir [82]. After a visual inspection of the year 2018 (Figure 33), the site showed a large meander bar (covering more than the extent of the original island) formed by the accumulation of sediments along the river, which was channelled to the north and east. Another limit of the single-pixel analysis is that it might return different values from the approach using polygons. In fact, the analysis of the site resurfaced extent is, of course, more accurate in determining which site is out of the water, because a single pixel might indicate a submerged site when only a small part of it has emerged. However, the zonal histogram analysis is limited by the difficulties in determining the extent of archaeological sites both from remote sensing and on the ground (Section 2.2), and it can reduce the samples significantly. Both analyses are also influenced by the fact that the images used are monthly or annual composites, and many variations might be observed with shorter temporal intervals, as shown by the 2018 monthly analysis.
Of course, any accuracy in measurements is tied to the precision at which it is possible to georeference archaeological sites in GIS. While extreme care was taken during this process (see Section 2.2) by refining the dataset and including polygons only with specific criteria, the impact of any inaccuracy should be acknowledged. This will manifest mainly in the pixel analysis, because if the actual position of a site is different, the result will change as well. Of course, accurate field data coupled with a thorough analysis of remote sensing or cartographic data, will improve any inherent inaccuracy in the data.
Another element to consider is the possible atmospheric impact on the results. All the images used to generate the NDWI were acquired at the TOA reflectance, meaning that they have not been atmospherically corrected. It must be said that the use of ratios (e.g., NDWI) can already lower the atmospheric impact [83][84][85][86] (see also Lillesand et al., 2015 [79], pp. 518-522). Applying a low cloud filter to generate cloudless images and a cloud masking function also minimises these effects [74,87]. Nonetheless, its impact should not be underestimated, because the two products will have a slight difference in the pixel values (atmospheric correction usually lowers reflectance [88][89][90]. It should be mentioned that, if applied to the identification of archaeological remains, BOA reflectance data are, however, highly suggested [91][92][93]. Future studies focused on single areas might use atmospherically corrected data and compare TOA data results to see which one works best, especially when reference ground data are not available.
Lastly, although it could be argued that the medium resolution of the images used could have a negative impact over the classification, the results highlighted above show that Sentinel-2-and Landsat-generated NDWI can still be used with a high degree of reliability. Of course, accuracy will always be tied to the resolution of reclassified and reference images, and higher-resolution data could always improve the results. However, as it has been already shown, when these are not available, testing with different thresholds can strengthen the results from medium-resolution images and eliminate noise from other ground sources caused by the surrounding context of each area [51,67].

Conclusions
Spectral indexes applied to medium resolution satellite images over a wide time frame can help in monitoring the effects of the receding waters on archaeological sites. The NDWI can distinguish between water and land surfaces in a useful way, allowing us to understand when a specific area is out the water or not. Applying the index on different images over time and by quantifying the results of the pixel analysis, it is possible to highlight years with very low or high water levels, with results partially comparable with the Global Surface Water (GSW) dataset [94] (Figure 3e), but also sites that resurfaced more often than others during the considered time frame. Knowing which sites are more likely to emerge can prove useful, because it is likely that these sites will resurface again if there is a decrease in water level in the reservoirs. This, in turn, can help in planning future targeted investigations. In fact, by knowing in advance which sites are more likely to emerge from the ebbing waters, salvage surveys, targeted excavations, or assessments of the site's conditions can have a more precise lead. On the other hand, knowing the number of emersions and coupling this information with the site status before its flooding (e.g., the type of archaeological evidence recorded) will help in a preliminary assessment of the site status before any future investigation.
The present methods and results, together with the difference between the single pixels and the zonal histogram analysis, also highlight the importance of accurate field recording of the extent of archaeological materials, and the need for more accurate and available geospatial data for archaeological sites. Moreover, field campaigns will help not only the cultural heritage management, but also the validation of the results.
In the end, the methodology presented here has the advantage of being easy to reproduce and apply to different areas, but also easily accessible by being based on free software and easy to access data, if we exclude archaeological sites data. In fact, the NDWI reclassification and analysis is faster and easier to reproduce than more in-depth methodologies such as the GSW, while still yielding promising results. Moreover, this approach can set a starting point to address a hitherto missing element of many salvage projects, i.e., the evaluation of the post-flooding impact of reservoirs on cultural heritage. In fact, apart from a few examples that inspired the development of this methodology [22,23], submerged sites have been considered mostly lost. However, it is clear that, at least for the lakes presented in this study, this may not always be the case.
Lastly, it has been already highlighted how costs, availability, and ease of access to commercial satellite images for the Near East are still a major hindrance for any research aimed at monitoring the archaeological cultural heritage, and how archaeologists had to cope with these restrictions in different ways (see, e.g., Rayne et al., 2020 [53] and related bibliography). In this sense, the daily availability of Sentinel-2 images and the possibility to monitor archaeological sites frequently, prove that these satellites images are and will be among the most useful free tools in any future remote sensing application.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2072-4 292/13/4/786/s1, Figure S1: Temporal coverage of satellite images used in the analysis; Figure S2: Emersion rate of archaeological sites in the Haditha Dam; Figure S3: Emersion rate of archaeological sites in the Mosul Dam; Figure S4: Emersion rate of archaeological sites in the Hamrin Dam. Code S1: Google Earth Engine code for the generation of the NDWI composite, QGIS models and scripts, and R code for images generation, pixel and zonal histogram analysis are available on Github: https:// github.com/andreatitolo/IraqEmerginSites, and archived on Zenodo (doi:10.5281/zenodo.4446664) NDWI composites, due to their size, are available on figshare (links for the annual and monthly NDWI composites are provided in the Github repository). Code S2: Link to Google Earth Engine repository: https://code.earthengine.google.com/?accept_repo=users/sapienza_at/IraqEmergingSites.

Data Availability Statement:
The data presented in this study are openly available on Github: https: //github.com/andreatitolo/IraqEmerginSites, and archived on Zenodo (doi:10.5281/zenodo.4446664) NDWI composites, due to their size, are available on figshare (links for the annual and monthly NDWI composites are provided in the Github repository).