Regional-Scale Systematic Mapping of Archaeological Mounds and Detection of Looting Using COSMO-SkyMed High Resolution DEM and Satellite Imagery

Abstract

“Tells” are archaeological mounds formed by deposition of large amounts of anthropogenic material and sediments over thousands of years and are the most important and prominent features in Near and Middle Eastern archaeological landscapes. In the last decade, archaeologists have exploited free-access global digital elevation model (DEM) datasets at medium resolution (i.e., up to 30 m) to map tells on a supra-regional scale and pinpoint tentative tell sites. Instead, the potential of satellite DEMs at higher resolution for this task was yet to be demonstrated. To this purpose, the 3 m resolution imaging capability allowed by the Italian Space Agency’s COSMO-SkyMed Synthetic Aperture Radar (SAR) constellation in StripMap HIMAGE mode was used in this study to generate DEM products of enhanced resolution to undertake, for the first time, a systematic mapping of tells and archaeological deposits. The demonstration is run at regional scale in the Governorate of Wasit in central Iraq, where the literature suggested a high density of sites, despite knowledge gaps about their location and spatial distribution. Accuracy assessment of the COSMO-SkyMed DEM is provided with respect to the most commonly used SRTM and ALOS World 3D DEMs. Owing to the 10 m posting and the consequent enhanced observation capability, the COSMO-SkyMed DEM proves capable to detect both well preserved and levelled or disturbed tells, standing out for more than 4 m from the surrounding landscape. Through the integration with CORONA KH-4B tiles, 1950s Soviet maps and recent Sentinel-2 multispectral images, the expert-led visual identification and manual mapping in the GIS environment led to localization of tens of sites that were not previously mapped, alongside the computation of a figure as up-to-date as February 2019 of the survived tells, with those affected by looting. Finally, this evidence is used to recognize hot-spot areas of potential concern for the conservation of tells. To this purpose, we upgraded the spatial resolution of the observations up to 1 m by using the Enhanced Spotlight mode to collect a bespoke time series. The change detection tests undertaken on selected clusters of disturbed tells prove how a dedicated monitoring activity may allow a regular observation of the impacts due to anthropogenic disturbance (e.g., road and canal constructions or ploughing).

1. Introduction

Digital elevation model (DEM) is currently a common dataset in archaeological research. Archaeologists and remote sensing analysts use DEMs for basic display of the morphologic and environmental setting of the study area, as well as more elaborated visualizations of landscapes through derived products (e.g., hill-shading, slope analysis, and local relief modeling [1]). The latter are typically utilized to support geographic information system (GIS)-based studies aiming at interpretation of detected features and sites in their topographic context (e.g., [2,3]), predictive modeling (e.g., [4]), viewshed, least-cost path analysis (e.g., [5,6]), or mapping and analysis of ancient monuments, urban, hydraulic, and agricultural networks (e.g., [7,8]).

DEMs are most frequently generated from: (1) airborne sensors, (2) drones, (3) terrestrial laser scanners, and (4) satellite remote sensing (an overview of the semantic categorization of potentials and practicality of such platforms for elevation data acquisition is provided in [9]). Light detection and ranging (LiDAR) data collected from airborne sensors allows for wide area coverage [1]. Terrestrial laser scanners are instead used to collect bespoke data for specific archaeological projects [2], mostly for site recording before or during excavation or conservation (e.g., [10]). More recently, LiDAR scanning from drones has emerged as a new solution for topographic data capture in archaeology; however, certain challenges (e.g., improvement of detection success, accuracy, and cost-effectiveness) are yet to be fully addressed [11].

In parallel, since the early 2000s, space agencies have increasingly released global DEMs that were produced using the following main methods: (1) stereo-photogrammetry and (2) interferometric synthetic aperture radar ((InSAR) e.g., [12]).

In the first case, image data of the same area are collected from slightly different positions, with different viewing angles along the satellite orbit; thus, differences in elevation result in a different parallax that can be measured and converted into elevation differences [9]. Examples of such global DEMs are the Advanced Land Observing Satellite (ALOS) World 3D–30 m (AW3D30 [13,14]) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global DEM (ASTER GDEM, at 30 m resolution [15]). In particular, AW3D30 is generated from stereoscopic optical images that are processed by means of a software encompassing, in order, tie-point generation, image orientation, image matching for height-calculation, and masking of outlier areas. The full details are described in [13]. Similar workflow is used by Biswas and Pandy [16] to generate a DEM from the Indian Remote Sensing Satellite, known as Cartosat-1 (or IRS P5).

In the second case, the input data are two SAR images collected by active microwave spaceborne sensors, ideally at the same time (in practice, with the shortest time separation possible), from two different points along the platform orbit, exploiting the spatial separation (the so-called baseline) either between two SAR antennas or two vantage points of the same SAR antenna [17]. An interferogram is formed by cross-multiplying the reference SAR image with the complex conjugate of the second. Assuming that the imaged ground is not changed during the time separating the two SAR acquisitions [9], the interferometric phase is proportional to the signal travel path distance and hence can be used to derive the landscape topography [17].

SAR-based DEMs include the Shuttle Radar Topography Mission, SRTM, at 30 and 90 m resolution [18]. SRTM is the most renowned example of the so-called single-pass InSAR system, by which both SAR images are acquired at the same time. On the contrary, repeat-pass constellations, such as the Italian Space Agency (ASI)’s COnstellation of small Satellites for Mediterranean basin Observation (COSMO-SkyMed), which is the main data source for the present work, provide SAR image pairs collected by different satellites observing the same location with the minimum time-lag equal to one day.

Given that SAR sensors collect data regardless of weather and day or night conditions, DEM generation from InSAR is essential for applications in remote areas where the photogrammetric approach to DEM generation using optical imagery may be affected by unfavorable weather conditions, e.g., cloud coverage [17]. For further details, the reader should refer to the overview of the technical properties of these DEMs provided in [19].

The above cited ASTER, ALOS, and SRTM global DEMs proved advantageous, owing to free accessibility, global spatial coverage, and good height and horizontal accuracy. As soon as archaeologists became acquainted with these DEMs, a substantial body of literature was published to prove their benefits for archaeological and paleo-landscape research. For instance, Wiseman and El-Baz [20] illustrated a rich portfolio of examples relying on the 90 m resolution SRTM DEM. More recently, analysts have also developed DEM-based methods for improving the visual interpretation of landforms to overcome limitations of traditional photograph-based approaches and ease the implementation in GIS and remote sensing software (e.g., [21]).

As archaeologists work at a variety of spatial scales, i.e., from regional down to single excavated unit, global DEMs at 30 or 90 m resolution often lack the spatial resolution and accuracy needed for site-specific work [22]. Furthermore, at local scale the random presence of artefacts, data gaps, or processing errors may occasionally hamper the suitability of global products for archaeological research purposes, despite their nominal performance. Therefore, archaeologists are increasingly exploiting DEMs at enhanced resolution, purposely processed from satellite acquisitions tasked over specific areas of interest.

For example, Malinverni et al. [23] exploited stereo-pair processing of panchromatic and multispectral Pleiadés acquisitions with resolution of 0.5 and 2 m, respectively, to produce the first detailed topographic documentation of Khinis in Iraqi Kurdistan. Gupta et al. [24], instead, used a 5 m grid spacing DEM generated from Cartosat-1 stereo pairs to locate archaeological mounds in the surroundings of the Buddhist site of Sannati in India and traced their elevation profiles.

In the SAR domain, a technological revolution was introduced by the TerraSAR-X add-on for Digital Elevation Measurements (TanDEM-X) satellite from the German Aerospace Center (DLR). Its primary objective to generate a global and high resolution DEM with unprecedented accuracy was accomplished in 2016 over 97% of the land masses, with 12 m ground pixel spacing at the Equator [25]. DEMs of higher resolution can be obtained from TerraSAR-X StripMap and High Resolution SpotLight imagery at 3 and 2 m, collected in alternating bistatic and monostatic acquisition modes. Erasmi et al. [26] used such DEMs to enhance the micro-topography of the ancient city of Magarsos in Cilicia, Turkey, while Rutishauser et al. [27] detected paleo-channels and assessed their impact on changes in the settlement patterns. The StripMap bistatic DEM also proved capable to improve the delineation of the eastern walls of the ancient city of Apamea in Syria, in the context of the local topography and modern features of the three dams confining the namesake lake [28].

A further space asset that can currently provide DEMs of potential use in archaeology is the COSMO-SkyMed, funded by ASI and the Italian Ministry of Defence as a dual-use (i.e., civilian and defense) end-to-end Earth observation system consisting of four X-band SAR satellites (i.e., CSK1, CSK2, CSK3, and CSK4) [29]. Their derived DEM products are generated by means of InSAR processing of co-registered scenes, in any acquisition mode (except for polarimetric ones), i.e., StripMap HIMAGE (SM) and Enhanced Spotlight, to cite the highest resolution ones that are fit for archaeological research purposes [30]. When atmospheric artefacts in the InSAR phase are properly mitigated, SM pairs with a 1-day temporal baseline and 150–200 m perpendicular baseline can lead to DEM less than 5 m absolute vertical accuracy [31]. The same SM imaging mode proved successful to create DEMs through radargrammetry, providing vertical accuracy comparable with that of TerraSAR-X StripMap DEMs [32].

Despite these achievements and the growing request of COSMO-SkyMed imagery from users of the cultural heritage sector [33], in the scientific literature, there are still few studies that exploit COSMO-SkyMed DEMs for archaeological research. Among these are the work by Salvini et al. [34], where COSMO-SkyMed SM DEMs processed with InSAR and radargrammetry allowed the authors to delineate a paleo-drainage in a wide area surrounding the Siwa and Al-Jaghbub oases of the western Sahara Desert between Egypt and Libya, as opposed to unsuccessful attempts relying on the SRTM and ASTER DEMs, and via photo-interpretation of high-resolution optical imagery.

This literature refers to the capability of COSMO-SkyMed as a SAR system to reveal buried or partially visible features. We are not aware of publications or studies demonstrating that COSMO-SkyMed images or derived DEMs have been applied by archaeologists to systematically record and map exposed archaeological deposits and features raising from the ground level, such as archaeological mounds.

The latter represent one of the fundamental categories of archaeological research [35], especially in the field of Near and Middle Eastern archaeology. Commonly referred to as “tells” by archaeologists, archaeological mounds are formed by human and geological activities and deposition of large amounts of anthropogenic material and sediments over thousands of years, after alternating phases of mud-based building construction, collapse due to abandonment, weathering or intentional destruction, reconstruction, and repairs. As such, tells are the most important and prominent features in Near and Middle Eastern archaeological landscapes [36]. The study of their spatial occurrence can reveal insights into the emergence, development, and organization of the first complex human societies [35].

Sherratt [37] was among the first scholars to identify the potential usefulness of satellite DEMs to search for archaeological tells, in the perspective of overcoming the known challenge of discovering where tells are on the ground in an accurate and precise way, despite consolidated robust methodologies for mapping and field-survey. Soon after, the pioneering study by Menze et al. [35] demonstrated that the 90 m SRTM DEM could help to scan areas that were not archaeologically surveyed on a supra-regional scale, and to pinpoint likely tell sites. In their test run with a machine learning classifier that identified probable tell sites from the terrain model, the authors found a sharp increase in the detection probability for sites with a height of more than 5-6 m, i.e., a value approaching the SRTM data accuracy.

The assumption that DEM of better spatial resolution and posting can allow a better visualization of mounds finds contrasting evidence across the literature, and its validity may also depend on data (post-)processing. For example, Gupta et al. [24] found that features with diameters approximately greater than 100 m can be visualized in the DEM, while smaller mounds were visualized in the 3D anaglyphs produced from the input images. Furthermore, spatial resolution is a trade-off with other parameters, e.g., spatial coverage. While mapping tells, archaeologists are required to maintain coverage of their survey, extensive enough to discern significant spatial patterning [38].

It is therefore with specific regard to the state-of-the-art and user requirements summarized above that this paper aims to demonstrate that DEMs generated from COSMO-SkyMed data can be a valuable dataset to conduct a digital survey of ancient mounds and features over wide remote geographic areas, under search and classification criteria complying with well-established methods of archaeological research.

The demonstration is undertaken across the Governorate of Wasit in central Iraq, where archaeological surveys have been undertaken in the past years, and current published compilations of site positions highlight significant knowledge gaps about the spatial distribution of tells and archaeological deposits. The study area also allowed testing the methodology of tell detection in conditions of limited to absent ground-truth validation, which at present is not an unusual situation in the Near and Middle East due to logistic difficulties, restrictions, and security concerns.

Accuracy of COSMO-SkyMed products is assessed against openly distributed DEMs (i.e., 30 m SRTM DEM and AW3D30), which are among the most commonly exploited by archaeologists worldwide. The results of the systematic mapping are then discussed in relation to the improvements produced in the knowledge of type, distribution, and condition of local archaeological deposits. This is crucial information given that archaeological heritage in Wasit is currently at risk of vanishing due to either natural erosion and weathering, encroachment of anthropogenic activities (e.g., ploughing, infrastructure, modern settlement, and dam construction), or looting. Through integration with other remote sensing data, COSMO-SkyMed products proved helpful to identify looting incidents and define a data analysis workflow that, in the future, may inform the compilation of a list of sites requiring priority attention by local stakeholders and authorities for heritage conservation.

2. Materials and Methods

2.1. Archaeological Setting of Wasit, Iraq

Figure 1 shows the location of the study area, in central Iraq, and the spatial extent of the bespoke collection of COSMO-SkyMed SM acquisitions that was used in this research (see Section 2.2.2). Using the Delmej reservoir as a key landmark, the satellite data were tasked to cover, on one side, nearly 7200 km² of alluvial plain from Al Wihda, 40 km south of Baghdad, up to Delmej, and on the other side, nearly 5500 km² of agricultural fields, from the reservoir to the east toward the Tigris river.

From an archaeological point of view, for millennia, Wasit has witnessed the rise and decline of different civilizations (e.g., Mesopotamian and Arabic). Use, abandonment, and reconstruction of settlements across time led to the formation of tells. The exploitation of the natural resources of the region left traces that can still be found, such as the so-called “hollow ways” formed by erosion from the traction of people and/or their animals walking along the same route repeatedly [39].

On-the-ground documentation of tells and archaeological deposits in Wasit (and the Qadisiyah to the south) dates back to the 1920s and Robert McCormick Adams and colleagues’ surveys [40] that were undertaken in the 1960s and early 1970s [41]. This was also the period when a long-term project of economic development was undertaken in Iraq and the ~616 km² Delmej reservoir was created to contrast the risk of desertification and provide a source of water to the local community. This infrastructure significantly impacted the conservation of cultural heritage, as demonstrated by a recent multi-temporal satellite imagery-based analysis published by Marchetti et al. [42] and Zaina [43], in the framework of the QADIS survey project investigating the south and eastern region of Qadisiyah as a joint Iraqi-Italian initiative of the Alma Mater Studiorum-University of Bologna and the State Board of Antiquities and Heritage (SBAH) of Iraq [44]. At the Delmej reservoir, Jotheri et al. [39] also documented hollow ways, of which some are submerged under the shallow water of the current marshland, and others are found at the edges of the marsh.

While these studies confirm that the archaeological community is active in unveiling the full archaeological potential of this region and reconstructing the organization of the first complex human societies, on the other side, they highlight the existing knowledge gaps across wide strips of Wasit (see also Section 2.4). Given that the entire region has not yet been entirely systematically surveyed with satellite imagery and derived products, further threats to conservation (and thus toward the compilation of a comprehensive and accurate listing of sites) include looting [45], agricultural activities [42], and economic development-related activities, alongside natural processes [43]. Agricultural and construction activities, in particular, are considered an under-analyzed phenomenon that continues spreading in absence of an effective control by the authorities [42].

2.2. Digital Elevation Model (DEM) Data

2.2.1. Global Medium Resolution Data

To cover the entire study area, NASA’s 30 m resolution SRTM DEM ([18]) was exploited. In particular, we used the SRTM GL1 Global 30 m Ellipsoidal (SRTM_GL1_Ellip) DEM distributed by OpenTopography, i.e., the product that was obtained by conversion from the orthometric version using the EGM96 geoid model [46]. While its mission requirement was to achieve absolute height accuracy below 16 m, validation showed that in Eurasia the 90% absolute height error is of 6.2 m [47]. From the same repository, the ALOS World 3D-30 m Ellipsoidal (AW3D30_E) was accessed [48], i.e., the product obtained by the same conversion from the original orthometric version [13,14]. The latter is characterized by height accuracy of 5 m [14]. The results of the comparison between SRTM_GL1_Ellip, AW3D30_E, and COSMO-SkyMed SM DEM to assess the accuracy of the latter are presented in Section 3.1.

Additionally, the TanDEM-X 90 m DEM was retrieved. This product is derived from the global DEM with a 0.4 arcsec (12 m) posting and has a reduced pixel spacing of 3 arcsec, which corresponds to approximately 90 m at the equator [25,49]. Given such a low resolution compared to the other DEMs, the TanDEM-X 90 m DEM was only used for a regional visual comparison of topographic features but not for the accuracy assessment. In this regard, better TanDEM-X products are provided by DLR (e.g., the original 12 m version), but their provision is under the charge of a service fee and therefore is outside the remit of this research.

For sake of completeness, we finally mention that the 25 m resolution SRTM X-band DEM by DLR [50] would have been another ideal product to compare with COSMO-SkyMed SM DEM, given that they are both derived from X-band SAR data. However, it covers only a marginal western portion of the study area, therefore insufficient for a meaningful comparison.

2.2.2. High Resolution COSMO-SkyMed Data

StripMap HIMAGE HH-polarized images, with 3 m resolution (single look) and 40 km swath width, were collected by the SAR sensors onboard COSMO-SkyMed satellites in right looking mode and with 3.1 cm wavelength (i.e., 9.6 GHz frequency) and were used as input for the generation of DEM products for the study region.

The image processing for DEM generation was conducted based on InSAR (e.g., [12]). This technique has often been used for high-resolution topography estimation from pairs of SAR images (e.g., [51,52,53]), including COSMO-SkyMed data (e.g., [34,54]).

When identifying the SAR pairs to generate the DEM for Wasit, two main technical aspects were accounted for. First, short temporal baselines B_T were required to minimize the potential effect of phase noise due to temporal decorrelation (e.g., [55]), especially across the cultivated fields spread across most of the study region. Therefore, the COSMO-SkyMed constellation was purposely tasked to acquire tandem pairs, with B_T of 1 day. This configuration addresses the technical requirement of the shortest B_T to minimize the likelihood of coherence loss, i.e., a phenomenon that can interfere with the generation of a data gap-free DEM. Until May 2019, the tandem configuration was achieved by CSK2 and CSK3 satellites, while starting from June 2019 by CSK4 and CSK2 [30].

Secondly, only large perpendicular baselines B_P could allow for precise estimation of the topography, given the inversely proportional relationship between the standard deviation of the InSAR-derived height estimates and B_P (e.g., [56]). However, an upper limit at the critical perpendicular baseline B_Pc was also set to avoid full loss of coherence due to spatial decorrelation.

Key information on the tandem pairs used for Wasit is summarized in

Table 1. COSMO-SkyMed StripMap HIMAGE scenes used for the InSAR processing and generation of the five DEM products. Notation: asc.: ascending; desc.: descending; θ: look angle at scene center; β: incidence angle at scene center; B_P: perpendicular baseline; B_Pc: critical perpendicular baseline for flat terrain; h_a: height ambiguity; #: tile ID number.

Tile ID	Input Dates	Beam	θ [°]	β [°]	Orbit	BP [m]	BPc [m]	ha [m]
#1	17/06/2018–18/06/2018	H4_04	29.0	32.2	asc.	−856	2027	6.2
#2	16/05/2018–17/05/2018	H4_03	26.5	29.3	asc.	+976	1782	4.9
#3	25/04/2018–26/04/2018	H4_07	33.6	37.4	asc.	−1043	2551	6.1
#4	03/02/2019–04/02/2019	H4_0A	18.1	19.9	desc.	−668	1100	4.7
#5

. Three pairs were acquired along ascending orbits (at ~05:40 local time) to cover the central and eastern sectors of the study region, i.e., including the Delmej reservoir and the plain between al-Kut and Al-Hay. Another pair of scenes was acquired along descending orbits (at ~18:00 local time) to cover the wide alluvial plain west and northwest of the Delmej reservoir. While the former scenes were acquired with a 50 km length in azimuth, the latter were sensed with a total length of 78 km to allow an enhanced coverage of the study area, and later cut into two 45 km-long frames.

Although ground-truth validation could not be collected and thus made available for this research for obvious reasons of accessibility constraints due to security issues and, later, due to the COVID-19 pandemic, the COSMO-SkyMed constellation was tasked to collect acquisitions at the same time as the originally planned field surveys. This is a change of paradigm, given that published studies using satellite-based DEM products for archaeological survey mainly relied on DEMs that were collected years before the actual archaeological survey. It is reasonable to assume that the major topographic features do not change in such a short time interval, and therefore older DEMs can still be suitable for an archaeological survey. However, when the focus is to capture archaeological features at a finer resolution and scale in a dynamic and highly changeable landscape, this temporal aspect may play a relevant role. As such, tasking satellites such as the SAR that are able to operate under basically any weather conditions allows for a more suitable timing of data capture that, in archaeology, is commonly achieved by means of airborne LiDAR or drones (though with generally much higher costs).

The DEM products after InSAR processing were generated in GeoTIFF format (float, 32 bits per sample), and projected to the WGS84/UTM zone 38N coordinate system. Their posting (pixel spacing) is 10 m by 10 m, and their vertical datum is the WGS84 ellipsoid. Relative and absolute height and horizontal accuracies of each DEM tile meet (or exceed) the specifications as per the COSMO-SkyMed mission and products description document [57]. The latter indicates a relative height accuracy ≤6.9 m for DEMs generated using SM SAR pairs with B_P between 100 and 3500 m (with an accuracy of 20%) and coherence of at least 0.8. While these figures provide an indication of the achieved quality standard as per the product handbook, the specific values of accuracy of each product depend on the characteristics of the SAR pair used for its generation and the coherence of the scene.

The so-called “height ambiguity” (h_a) of each SAR pair, for instance, provides an indication of its interferometric performance in estimating the topography of the observed region (e.g., [56]). This parameter indicates the altitude difference, leading to a 2π interferometric phase change (after flat-Earth removal), and can be computed as follows:

$$h_a=\frac{\lambda H\tan \theta }{2B_P}$$

(1)

where H is the satellite orbit altitude (i.e., 619 km for the COSMO-SkyMed constellation), and B_P is the perpendicular baseline of each processed pair. As a general rule, the higher the B_P, the smaller the h_a, and therefore the more accurate the height measurement is. As shown in Table 1, the data pairs used for the analysis allowed for h_a of 4.7 to 6.2 m, which indicates a very good sensitivity to topography for the entire set of DEM products, in compliance with the nominal height accuracy standards specified in the product handbooks.

In order to estimate the accuracy of the InSAR-derived DEMs, the elevation values from the COSMO-SkyMed DEM tiles (h_CSK) at each location i were compared with those retrieved at the same location within the 30 m spatial resolution SRTM DEM (h_SRTM) and the ALOS World 3D—30 m DEM (h_ALOS) products. The elevation differences Δh_i were then calculated as follows:

$$\mathsf{\Delta }{h}_i=h_{CS{K}_i}-h_{SRT{M}_i}$$

(2)

$$\mathsf{\Delta }{h}_i=h_{CS{K}_i}-h_{ALO{S}_i}$$

(3)

Quantitative comparison and statistics calculations to assess COSMO-SkyMed SM DEM accuracy were made with regard to the AW3D30 and SRTM at 30 m resolution products, as these are the typical DEMs used in the literature by archaeologists investigating Asian arid environments and landscapes and, as such, were preferred versus the 30 m resolution ASTER GDEM. The TanDEM-X DEM at 90 m resolution was also not considered for quantitative comparison because of its much lower spatial resolution.

2.3. Other Satellite Imagery

To support the recognition of archaeological features, we used Google Earth and ESRI World Imagery, the latter accessed through ArcGIS software and dated on 2015. However, low or degraded spatial resolution, as well as the lack of time lapses, partially hampered the visibility of features in Google Earth across most of the study area.

Since agricultural and anthropogenic activities have significantly impacted the landscape (see Section 2.1), we used declassified CORONA satellite imagery from the United States intelligence satellite program [58]. The CORONA mission was launched from Vandenberg Air Force Base in 1960 and operated until May 1972. Its image collections were declassified in 1995 [59]. Of the various CORONA systems, the images used in this research were acquired from KH-4B in 1968 and are identified with unique image designators from ds1107-2170da122 to ds1107-2170da132. In particular, these images consist of black and white photographs collected through two panoramic cameras with a separation angle of 30°, with one camera looking forward and the other looking aft, on long film strips, each covering 8.6 by 117 km swath, with a nadir resolution of approximately 6 feet (i.e., ~1.8 m) [60]. Such high resolution was suitable to enable us to document land use and archaeological deposits at the time of image collection (i.e., 1968) and compare with the current visibility and detectability in the COSMO-SkyMed DEM.

Figure 1 and Figure 2a show the spatial extent and coverage of the tiles that were accessed as orthorectified GeoTIFF files (native Spatial Reference System: EPSG 4326/WGS84) from the CORONA Atlas of the Middle East [61]. We sourced these files as already orthorectified products. Full details about the methods by which CORONA images were efficiently and accurately orthorectified and, as such, are distributed through the online platform, are reported in [60]. According to Casana et al. [62], the accuracy of most of the orthorectified images ranges from 3–10 m at nadir to 20–80 m at the edges. Such accuracy was suitable to compare with the 10 m horizontal posting achieved by the COSMO-SkyMed DEM products. Visual inspection of the CORONA images at the local scale (taking precise geographic locations and toponyms that have not changed over time as spatial references) over the study area did not highlight geometric distortions or geolocation errors that could affect the following analysis and interpretation.

To assess the present-day conditions of the mapped archaeological features, we also referred to the Copernicus Sentinel-2 multispectral image collection available over the study area in the period 2017–2020, i.e., coeval to the COSMO-SkyMed SM acquisitions to generate the DEMs and the planned field activities. According to the methodological approach described in [63,64], out of the 13 bands provided by Sentinel-2 data we mainly focused on the visible and near-infrared (NIR) bands at 10 m spatial resolution and combined them to generate false-colored infrared composites and basic spectral indices (e.g., Normalized Difference Vegetation Index–NDVI). Aligning with methodologies published by [38,65], Sentinel-2 imagery was also helpful to compare the spectral profiles of the mapped features with their radar backscatter signal in the COSMO-SkyMed SM data, and thus corroborate their identification and type classification.

2.4. Ancillary Data

As well demonstrated in the literature [69], a reference source for systematic archaeological survey of earthen mounds in Asia is provided by the Soviet military maps of the 1950s [66,67], owing to their precision and the early date at which they were produced [69]. To correctly interpret the topographic anomalies drawn in these maps and check the nature and type of archaeological mounds (Figure 2c), we referred to the legends translated in English [68] and in particular to the symbols of “62. Burial mound (height indicated in meters)”, “456. Triangulation point on burial mound”, “472. Burial mound”, and “473. Tailings pile” (Figure 2d). In several instances, the maps not only report the toponyms, but also the explicit indication of “Телль”, i.e., “tell” in Russian (Figure 2c). The term is often preceded by the Russian term “бугoр”, i.e., “bugor” (“hillock”; abbreviation: буƨ; Figure 2c). In many other cases, these indications are absent, but this could be expected given that mapping of archaeological mounds was not one of the goals of such cartographic efforts [69]. However, surveyors differentiated artificial from natural relieves by drawing the former with hachured contour lines (pointing toward lower elevation), and the latter with dashed lines. This is helpful information that we have accounted for (although requiring verification), in addition to the indication of the height value indicated at some of these tells. As demonstrated by Menze et al. [35], this value coincides with the height of the tell measured from base to top, and can be used for comparison with the tell height as observed in the DEM data. As with CORONA imagery (see Section 2.3), visual inspection of the Soviet maps at the local scale did not highlight any significant geometric distortions that may affect our assessment.

Figure 3 shows the spatial distribution of placemarks of Ancient Near East (ANE) sites that Pedersén [70] identified and collected in an open-access database, through visual identification of sites on Google Earth satellite images. The database has some declared limitations (e.g., towns smaller than 1-2 ha have often not been included, unless the site was already well-known or of special interest for a certain reason; some geographical areas and time periods have been more completely treated than others [70]). However, it provides an excellent entry point resource for the study of archaeological mounds in Wasit and suggests that the use of DEMs for purposes of regional surveying would be beneficial. While a good density of ANE placemarks is found over the study area stretching from south of Baghdad up to the Delmej reservoir, the much smaller number of mapped sites across the agricultural fields in the south-eastern part suggests a knowledge gap about the actual distribution of archaeological mounds.

Ancillary data on local hydrography (Figure 1 and Figure 3) include: waters of the Mesopotamian floodplain by Pedersén [71] representing modern waters systems (Euphrates and Tigris rivers, Drain), Islamic canals, ancient and geological waters and 1:50,000 scale Iraq-Water Courses (Rivers and Streams) by the United States National Imagery Mapping Agency (NIMA) [72]. Furthermore, we referred to Iraq waterways, alongside Iraq buildings, natural features, roads and places by OpenStreetMap.

2.5. Archaeological Mapping of Tells in Wasit

Topographic features and anomalies visible in the CORONA images and the COSMO-SkyMed DEMs were mapped and interpreted in a GIS environment. Given that the specific aim of the present work (see Section 1) was to stimulate further exploitation of COSMO-SkyMed DEM products in archaeology, we developed our use-case demonstration by considering the most probable and simplest of the exploitation scenarios, i.e., systematic visual identification and manual mapping of features under search and classification criteria complying with well-established methods of archaeological research.

CORONA, satellite imagery, and COSMO-SkyMed DEMs were overlaid by a 10 km × 10 km search grid, and each subset was systematically inspected on screen to map anthropogenic or natural features of interest within it, each with its correlated metadata. As far as the anthropogenic features, recorded ancillary data include annotations on the nature of the trace and an Index of archaeological reliability [73]. Each entry was assigned an ‘archaeological reliability’ score according to a five-tiered scale from 1 (very low) to 5 (very high). Reliability is defined here as the level of likelihood of the trace to have an anthropogenic origin. The score was assigned based on the number of repeated identifications over multiple examined sources. As a higher archaeological reliability would correspond to priority for the location being inspected during the ground-truth phase, this value was translated into a ‘Priority’ value. This approach enabled managing situations when the anthropogenic origin of a mapped object was uncertain or questionable, and all the data sources did not provide a specific indication (see also Section 2.4).

When allowed by the spatial resolution of the satellite imagery available, the mapping exercise also encompassed assessment and recording of looting traces.

3. Results

3.1. COSMO-SkyMed DEM Coverage and Accuracy

The five COSMO-SkyMed SM DEM products cover a total area of 9000 km², with elevations above the reference WGS84 ellipsoid ranging between −7.74 and +40.67 m (Figure 4 and Figure 5). The quantity of data gaps due to phase decorrelation span from 3.3% in tile #1 to 23.5% in tile #3 (

Table 2. Main characteristics of the five COSMO-SkyMed SM DEM products and accuracy statistics from their comparison with the SRTM_GL1 and AW3D30 DEMs at 30 m resolution. Notation: h_min and h_max: minimum and maximum elevation above reference ellipsoid; Min.: minimum; Max.: maximum; Av.: average; StD.: standard deviation; RMSE: root mean square error.

Tile ID	hmin [m]	hmax [m]	No Data [%]	SRTM_GL1 Δhi [m]				AW3D30 Δhi [m]
				Min.	Max.	Av. ± StD.	RMSE	Min.	Max.	Av. ± StD.	RMSE
#1	−7.74	+13.17	3.3	−26.02	+27.09	−0.04 ± 1.28	+1.28	−17.06	+50.76	−0.46 ± 1.41	+1.48
#2	−4.07	+21.12	10.0	−43.74	+40.29	−0.23 ± 2.64	+2.65	−17.69	+15.56	−0.13 ± 2.32	+2.32
#3	−1.52	+25.66	23.5	−29.42	+24.43	−0.39 ± 1.92	+1.96	−16.94	+25.83	−0.14 ± 1.39	+1.40
#4	+14.55	+40.67	12.0	−26.05	+13.48	−0.16 ± 1.69	+1.69	−12.64	+10.71	−0.04 ± 1.23	+1.23
#5	+15.42	+38.15	18.8	−22.66	+18.47	−0.16 ± 1.74	+1.75	−26.41	+42.27	−0.34 ± 1.42	+1.46

), reflecting the land cover and land use and, in particular, the Tigris river course, artificial channels across the large extents of cultivated fields (Figure 4b), and marshes and water bodies (e.g., Delmej marsh).

The observed minimum, maximum, average, and standard deviation of the elevation differences Δh_i for each of the five COSMO-SkyMed SM DEM tiles with respect to the two external DEMs (i.e., SRTM_GL1 and AW3D30) are summarized in Table 2, along with the observed root mean square errors (RMSE). The comparison is based on the number of valid DEM pixels (i.e., those including a height value) across the entire extent of each tile, namely ~19.3 million pixels for tile #1, ~18.3 for tile #2, ~15.3 for tile #3, ~16.6 for tile #4, and ~15.5 for tile #5.

Elevation values are well-aligned with those of the external DEMs, as the elevation differences are smaller than 0.5 m on average, with standard deviations typically below 2 m. The RMSE is also in the range of 1.2–2.6 m. For most tiles, the elevation differences appear slightly lower against the ALOS DEM, rather than the SRTM, though no systematic shift is observed through the tiles. These statistics provide a confirmation of the very good accuracy of the 10 m resolution InSAR-derived COSMO-SkyMed SM DEMs compared with the two global products at 30 m resolution, and their suitability to support the archaeological mapping of deposits and tells.

3.2. Visibility and Spatial Analysis of Archaeological Deposits

The different spatial resolution of the products (10 m vs. 30 m) implies an improved capability (by a factor of 9) of the high-resolution COSMO-SkyMed SM DEM to capture topography compared to the medium resolution external DEMs. In an earlier publication [30], we commented about the improved capability in the detection of archaeological tells that can be achieved with the 10 m posting of COSMO-SkyMed SM DEM through an example in the central agricultural landscape east of the Delmej reservoir. Similar improvement is achieved in various sectors across the wide study area, where temporal coherence between SAR image pairs was high and decorrelation due to intense agricultural activity did not cause data gaps in the final DEM.

The benefit toward better feature detectability can be more appreciated if the DEM product is compared not only with the input SAR scenes (displayed as geocoded and terrain-corrected amplitude images) from which the DEM has been generated, but also with the Copernicus Sentinel-2 multispectral images and the very high resolution (VHR) ESRI World Imagery basemap.

Figure 6 and Figure 7 show two examples of tells, the location and presence of which were known from field surveys. In the first case (Figure 6), the topographic pattern and associated profile extracted from COSMO-SkyMed SM DEM complement the clear signals in the visible, NIR, and radar bands that already suggested the presence of a mound with 4 m height and over 200 m diameter. These morphological properties match with those typical of archaeological tells [35].

On the contrary, in the second case (Figure 7), no evident surface mark or contrast in multispectral bands or distinctive radar backscatter versus the surrounding soil is found in satellite imagery. Therefore, the pattern identified in the COSMO-SkyMed DEM becomes essential to reveal the presence of the mound (Figure 7d). Additionally, the asymmetrical shape of the topographic profile (Figure 7e) suggests that the mound topography may have been partially flattened and altered over time. This also matches with the tire marks and signs of earth moving in the VHR basemap (Figure 7a), as well as the widespread smoothed reflectance in the Sentinel-2 false color composite (Figure 7b). Several other similar instances of archaeological tells difficult to detect in optical imagery were found across the entire study area.

Our findings achieved through data integration confirm what was observed in the specialist literature [35,37], i.e., that marks in multispectral imagery can be too unspecific to serve in detecting tells. Furthermore, anthropogenic activities and agricultural practice may have damaged or even leveled the mounds. Therefore, even if mounds were formerly characterized by specific marks, these may have been altered. In these circumstances, if there are no other specific environmental or soil conditions that can highlight surface anomalies relating to the presence of mounds [65], only optimally processed high resolution topographic products [35] can help to identify mounds.

COSMO-SkyMed SM DEM proved a valuable resource to achieve an overall picture of the spatial distribution of tells and archaeological deposits and start to make some considerations on their spatial patterns. In general, tells appear to be distributed along water courses that (apart from some exceptions) are mostly paleo-channels detected by comparing CORONA, Sentinel-2, and modern satellite imagery. The mapped paleo-hydrographic features exhibit a width ranging from 150 m to 1.5 km.

Figure 8 shows a very representative example located in the central portion of the north-western sector of the study area. The mounds as mapped through the Soviet maps (red polygons) systematically overlap with paleo-channels. This evidence suggests the existence of some relationships between the complex local hydrographic network and settlement patterns developed over time. Figure 8 also allows noting a discrete and widespread number of sites, the presence and location of which were verified and confirmed from the Soviet maps and CORONA imagery to COSMO-SkyMed SM DEM, albeit the different priority assigned during the mapping exercise.

Similarly, numerous correspondences between the analyzed datasets were found in the south-eastern sector of the study area (Figure 9), with an apparent higher density of sites at higher level of priority (i.e., 4 to 5) and larger diameter (Figure 10a and Figure 11a,b).

While a settlement pattern analysis and the archaeological interpretation of the mapped features are beyond the scope of the present study, from a statistical point of view, of all the mapped features, nearly 60% are classified from high to very high priority, while only 3% and 12% are of very low to low priority, thus of more probable natural origin and not archaeological features. The diameter of the mapped features is as wide as 100–200 m for 54% of the mapped tells, while it is less than 100 m for 15% of the sites (Figure 11b). Of the mapped features, 84% were equally visible in both the COSMO-SkyMed DEM and any of the other ancillary or satellite data (e.g., CORONA, Soviet maps, Sentinel-2 imagery), while 12% were detected in the COSMO-SkyMed DEM only. Furthermore, 22 entries (of which 13 and 4 are at very high and high priority, respectively) match with ANE placemarks, whereas the rest of the mapped entries are new and previously unmapped.

3.3. Condition Assessment and Impact of Looting

Google Earth and ESRI World Imagery combined with the location of the sites derived from COSMO-SkyMed SM DEM have clearly shown an impressive concentration of sites with looting evidence in the southeastern area of the study region (Al-Ahrar sub-district, administrative code IQG18Q03N01; Figure 10b). With the exception of the very few instances in which the evidence of looting is uncertain and requires inspection through more adequate spatial resolution imagery (Figure 11c), 42% and 38% of the identified looted sites are classified as very high to high priority, respectively (Figure 11d). The possibility of these hypothetical archaeological sites to genuinely be tells is demonstrated by the existence of looting traces in their proximity, as it is likely that the archaeological potential and historical value of the sites were known to local traffickers. To further corroborate this observation, of the 22 entries matching with ANE placemarks, 15 sites show clear evidence of looting, 3 require further verification, and only 4 show no evidence in the available data.

The second area most rich in looted sites is located toward the northwest (Markaz Al-Noamaniya sub-district, code IQG18Q03N02), but in this case the lower resolution of the satellite imagery partially hampered its visibility, therefore the records may be affected by inevitable under-estimations.

Figure 12 reports a selection of sites where, alongside signs of looting, VHR Google Earth and ESRI base map imagery and the color contrast allow the recognition of traces of structures, mostly in the form of walls and enclosures. In Figure 12a,b, looting pits appeared fresh on 1 January 2005 and faded on 13 March 2010. The more recent image better highlights the presence of anthropogenic structures and enclosures and suggests that some sites were not affected by looting. In Figure 12c,d, the image collected on 21 November 2004 shows the occurrence of looting activities and that they encroached a previously well-preserved structure.

4. Discussion

4.1. Assessment of COSMO-SkyMed DEM

The statistics calculated by comparison with SRTM_GL1 and AW3D30 DEM at 30 m resolution, with elevation differences being smaller than 0.5 m on average and standard deviations typically below 2 m (see Section 3.1, Table 2), reassure about the accuracy of the InSAR-derived COSMO-SkyMed DEMs. Alongside the improved spatial resolution owing to the 10 m posting (Figures 4b and Figures 5b), such properties prove COSMO-SkyMed SM DEMs as a valuable resource for archaeological mapping compared to the global DEMs. Further validation of the topographic height of the mapped tells may be accomplished via ground-truth survey data (e.g., ground control points surveying with kinematic GPS).

Our experiment demonstrates that it is possible to task satellite SAR acquisitions to be collected in parallel with planned field surveys, and thus generate a fine resolution and accurate DEM according to a functionality that (albeit with better spatial resolution, at higher costs) can be nowadays achieved with airborne and drone platforms only, but may be impractical and unfeasible in the Middle East and other inaccessible areas of the planet. In particular, the tandem COSMO-SkyMed SM pairs collected over Wasit met the technical requirements for DEM generation, i.e., the shortest possible temporal baseline of 1 day and acceptable range of perpendicular baselines (from 668 to 1043 m; see Section 2.2.2, Table 1). While the latter parameter contributes to reduce the height ambiguity (see Equation (1)), the former allows for lower probability of temporal decorrelation.

This aspect is of particular relevance for applications in a landscape such as that of Wasit to account for the possible interferences due to local land cover and anthropogenic actions. The high temporal dynamicity of agricultural practices and irrigation (with consequent changes of roughness and soil moisture content) make the landscape vary rapidly and this, in turn, causes temporal decorrelation that induces data gaps within the DEM. Therefore, the agility of COSMO-SkyMed constellation to be tasked to collect pairs and repeat acquisitions in short time periods is advantageous and, in the present research, was an essential feature of the constellation to decrease the amount of data gaps due to phase decorrelation up to 3.3% in a highly cultivated landscape in the easternmost sector of the study area. Marshland and water bodies will always be detrimental to DEM generation, as proved by the higher proportions of data gaps found within the swath of the SAR pairs including the Delmej reservoir (see Section 3.1; Table 2).

In summary, compared to the other space-borne DEMs, the use of COSMO-SkyMed images for DEM generation is beneficial because:

• At equal accuracy, the resulting DEM is of much smaller horizontal posting and therefore enhanced resolution;
• If seasonal and contextual anthropogenic factors are accounted for, acquisitions can be tasked to achieve a DEM representing the situation at a given time window (which may also match with the timings of the terrestrial surveys undertaken by field archaeologists, if allowed);
• New acquisitions can be tasked over time to generate new DEM products for topographic change detection analysis (although this was beyond the scope of this research).

At equal requirement of processing efforts (which largely depend on the availability of adequate software and operator’s expertise), COSMO-SkyMed DEMs are more cost-effective than LiDAR and drone-based DEMs because the production cost of the input SAR images is much lower, up to no cost in cases of data exploitation for scientific purposes through announcements of opportunity and open calls for data access (see [30]).

4.2. Spatial Distribution of Archaeological Deposits and Tells

As explained above, tells and archaeological deposits were mapped across the study area of Wasit using COSMO-SkyMed DEM by means of a traditional methodology consisting of expert-led, intensive and systematic screening, visual identification, and manual digitization, and through the combination with a range of satellite imagery and historical records (see Section 2.3, Section 2.4 and Section 2.5). This methodological choice, alongside the fact that the exercise was intentionally run by an archaeologist instead of an image analyst, makes the present exercise be a realistic use-case proof-of-concept of what any other archaeologist or cultural heritage expert could develop if they used COSMO-SkyMed DEM products for similar study purposes, in comparable environments and archaeological landscapes. On the other side, future research may focus on testing automated or machine learning-based methods to speed up the mapping process. However, at this stage, such an advanced step would imply the following two conditions that are yet to be met: COSMO-SkyMed DEM products being already recognized and used as an established topographic dataset by archaeologists; computing skills being developed from the user side and adequate training samples being available for feature recognition and classification.

In the specific case of Wasit, COSMO-SkyMed DEM products proved advantageous not only to achieve a better visibility of topographic features and thus a better detectability (see Section 3.2), but also to overcome some drawbacks of the available global DEMs, such as artifacts affecting at regional scale the AW3D30 product. Compared to the ANE placemarks dataset, our results not only confirm the presence of previously mapped tells, but more importantly expand the tell inventory, especially in the southeastern region where more data gaps were present in the ANE database (see Figure 3 and Figure 9). Consequently, the knowledge about the spatial distribution of tells across the Wasit landscape has been significantly improved and some initial hypotheses about the association between the location of tells and local hydrographic networks can be formulated.

However, as correctly stated by Menze et al. [35], settlement and feature patterns are largely a reflection of the intensity of investigation, and such rules also apply to the present mapping exercise. Therefore, we cannot expect that this technological product (and its future analysis through automation and machine learning) can replace the traditional process of conducting surveys on the ground which, for a long time, has been the main approach for archaeologists to interpret the landscape topography and should continue to be as part of a complex physically engaging interpretive process [74]. On the contrary, we advocate for the use of bespoke DEM products such as those generated from COSMO-SkyMed SAR acquisitions to facilitate the detection and mapping of tells remotely, in order to better populate tell inventories and enrich the existing databases and, in turn, increase the samples to verify through on-the-ground surveys, when allowed, and to use for settlement pattern analysis and landscape topography interpretation.

It is with this scope in mind that the results are presented in this paper, alongside the other key aspect of “digital fieldwork”, i.e., documentation of archaeological features that are at risk of vanishing under the pressure of anthropogenic hazards.

4.3. Monitoring Strategies for Condition Assessment

The identification of sites that may have been leveled or the original morphology of which may have been altered, alongside the documentation of looted sites, are of outmost importance to understand the exposure of cultural heritage to hazards in a changing landscape, such as that of Wasit.

Our evidence of looted sites aligns with the evidence reported in the recent literature [42]. Furthermore, the spatial location of looting pits in proximity to or within yet unexcavated sites matches with on-the-ground observations that led archaeologists to infer that looters had a deep knowledge of the territory and knew where to search for goods [42]. Nevertheless, it is fair to acknowledge that the overall picture of spatial patterns gathered in this study may not be complete, given that some portions of the study area are covered by imagery at low or degraded spatial resolution, or that no imagery was recently collected, thus the observations were not up-to-date. Therefore, it is reasonable to imagine that more sites have been looted and the situation is further evolving in more recent times.

Furthermore, as recently reported in the neighboring area south of the Delmej reservoir [43], sites are also exposed to threats due to construction of roads, modern canals, and ploughing. Therefore, monitoring the site condition also means assessing the impact due to the anthropogenic use of the territory where the sites are located in.

A further step in this regard may be represented by monitoring selected areas that the mapping exercise has highlighted as at risk of being looting by means of tailored SAR acquisitions, collected on a regular basis at higher spatial resolution. Figure 13 shows an example of COSMO-SkyMed Enhanced Spotlight (SP) images at 1 m spatial resolution [30] that were tasked to be collected as a multi-temporal series over the period between 27 July 2019 and 24 March 2020, to monitor a 10 by 10 km wide portion of the landscape characterized by a high density of archaeological and paleo-environmental features and tells, for which the regional survey had highlighted several incidents of looting. Figure 13b highlights the observation scale that a COSMO-SkyMed SP image can allow and how it is possible to comparatively monitor different clusters of tells simultaneously.

Figure 14 shows the nearly 400 m wide Tulul el-Barakat TB1 site (id.57, known from [70]) that has been heavily looted before the end of 2004, as evidenced in Google Earth imagery (see red arrows in Figure 14a). The same looting pattern is still visible in the COSMO-SkyMed SP image (see red arrows in Figure 14b) thanks to the distinctive texture and appearance of looting pits, according to the SAR ascending observation geometry [75]. The COSMO-SkyMed SP image also highlights the areas excavated by SBAH in 2008–2010 [76]. Accounting for the precise spatial co-registration between the optical and SAR images, looting features appear not to have aged and vanished.

Figure 15, instead, provides evidence of four sites (id.66–69) that had not been already mapped before the present study and were found to be heavily disturbed by agricultural activities and narrow roads (see cyan arrows).

In both cases, we applied the change detection method according to the methodology developed for VHR X-band SAR imagery [75]. While the yellow-colored pixels mark areas of no change between 27 July 2019 and 24 March 2020, a substantial part of the Tulul el-Barakat TB1 site appears disturbed with patterns of decreased backscatter signals coinciding with either excavation or earth movement or ploughing (see, for example, the feature highlighted by cyan arrows in Figure 14c). Signs of a narrow pathway crossing the tell are visible (see black arrows in Figure 14c).

Similar evidence is found on the same date over sites id.66–69 (Figure 15b), with marked patterns along the canals and the narrow roads crossing or running near to the tells (see cyan arrows). However, some differences between the sites are found, with changes of higher magnitude in site id.66 versus 69.

It is well understood that this example is meant as a proof of concept only of a possible monitoring activity that may be used to inform local conservation practice, and hopefully translate into protection measures in the future.

5. Conclusions

The Governorate of Wasit in central Iraq encompasses an archaeological landscape that is still rich in tells and archaeological deposits that previous investigations have only partially unveiled. However, the current land use and anthropogenic activities pose serious threats to conservation, and therefore documentation of the surviving remains and archaeological features is required before they vanish. As such, this study region served as an ideal test-bed to experiment the COSMO-SkyMed radar imaging technology to generate high resolution DEMs and undertake a systematic regional-scale mapping exercise. Furthermore, given that tells in Wasit are frequently embedded among fields that are currently exploited for agriculture, the capability of retrieving accurate COSMO-SkyMed DEMs was also tested to counterbalance the known effects due to temporal decorrelation and consequent data gaps in the generated DEM.

Not only was the accuracy of COSMO-SkyMed SM DEMs quantified with respect to the SRTM and AW3D30, but also the tell detection and digital documentation exercise was undertaken through expert-led visual identification and manual mapping. In this manner, we emulated what a large percentage (if not the majority) of archaeologists would have likely done if they were provided with the COSMO-SkyMed SM DEM.

The demonstration highlights that, thanks to the 10 m posting and the consequent enhanced observation capability, the COSMO-SkyMed DEM is advantageous to detect both well preserved and leveled or disturbed tells, standing out for more than 4 m from the surrounding landscape. Through the integration with other optical products and historical maps, the COSMO-SkyMed DEM not only provides the confirmation of the spatial location of sites known from the literature, but also allows for an accurate localization of sites that had not been previously mapped.

On this basis, by analyzing the multispectral signature and radar backscatter signal from VHR optical imagery, Sentinel-2, and the COSMO-SkyMed radar amplitude images, we achieved a picture of the survived tells and of looting locations across the whole study area.

While we hope that the present study may stimulate further implementation of the COSMO-SkyMed technology and the associated mapping approach for similar systematic investigations in other geographic regions, we envisage that future developments may include: the experimentation of (semi-)automated methods for detection and mapping of tells (e.g., through machine learning, as tested by [77]), in situ validation and GPS surveying to assess the absolute vertical accuracy on height estimation (which in the present research was impeded by logistic and accessibility constraints, e.g., [78]), and a complete archaeological investigation and interpretation of the settlement patterns (e.g., [78], which is beyond the remit of this publication).

Finally, with regard to the possible use of COSMO-SkyMed DEMs for archaeological mapping in different cultural and geographic contexts, further progress may include testing the same type of SAR acquisitions, DEM generation, and systematic mapping in other locations and landscapes in Europe and beyond. In this case, we envisage that, alongside dedicated acquisitions, a search within the rich interferometric archives (collected also with 1-day temporal baseline) should be made. These image collections may provide observations with temporal and perpendicular baselines suitable to generate DEMs and document archaeological features in seasons with limited interference from agricultural and anthropogenic activities. These, indeed, represent potentially limiting factors for the generation of DEMs covering vegetated and cultivated study areas.