**1. Introduction**

The study area (Figures 1 and 2) lies in the fine-grained, low-gradient alluvial plain of the Piave River [1,2], which is part of the Piave megafan that has been forming between the Last Glacial Maximum and the Late Holocene [3]. Fluvial activity by the Piave River (a major Alpine river with a catchment in the Dolomites) and the Livenza River (a minor Prealpine river) ended here before the Roman Age [4], so in this time, the area corresponded to a stable alluvial plain. Due to relative sea-level rise in the Middle Ages, the area became part of an extensive system of coastal wetlands that was completely reclaimed only between the 19th and the 20th century. The current topographic surface, artificially drained, is 2 m below the sea level.

**Figure 1.** (**a**) Location of the study area. (**b**) Agea 2012 digital color orthophoto: (1) the study area; (2) the site investigated with geophysical prospections, core drilling, and stratigraphic analysis of the exposed sections.

In Roman times, the area of Ceggia was included in the southern part of the "*ager*" of *Opitergium* (currently Oderzo), not far from the inner edge of a sizeable lagoon-marsh extended in the vast plain between the ancient river courses of the Piave and Livenza (Figure 2a). The Roman road so-called *Via Annia* [5–7] represented the principal axis of communication on the eastern coast, marking the limit among different landscapes. The *Via Annia* crossed the territory connecting the Roman centers of *Altinum* and *Iulia Concordia* up to *Aquileia*. Different archaeological findings during the times in the area (Figure 2b) testify to its Roman occupation [8–10]. No details about the population or the layout of these landscapes are registered. The only information available is related to the so-called stream Piavon-Canalat that flows in this sector, and probably representing from ancient times one of the routes of communication and commerce [11,12] between *Opitergium*, the lagoon, and the Roman harbor to the sea, the so-called *Portus Liquentiae* (Figure 2a).

**Figure 2.** (**a**) Historical-topographical map; (**b**) Geoarchaeological map with traces from aerial photographs. Legend: (1) Roman cities; (2) Early Medieval towns; (3) probable site of the Roman port; (4) Roman roads; (5) main watercourses; (6) contour lines with equidistance 2 m; (7) southern countryside of *Opitergium*; (8) mechanical drainage areas; (9) extension of marshy-lagoon areas from 17th century cartography; (10) study area; (11) area studied with geophysical surveys, core drilling, and cleaning of the exposed sections; (12) paleochannels; (13) ancient lagoon canals; (14) loamy soils; (15) clayey soils; (16) archaeological findings known from the bibliography [8,10] and from data archive of the SABAP modified by the authors; (17) anthropic aerial photo lines; (18) sites of particular interest: A-Roman bridge; B-probable rustic settlement. The geomorphological elements (12), (13), (14), and (15) are taken from [1].

In the early Middle Ages, Ceggia was once again a marginal reality or rather a borderland, if we consider that it must have been on the boundary between the Byzantine territories of the coast and the Longobard hinterland, as well as between the mainland and spaces dominated by marshy and lagoon waters. Di fferent studies conducted in the area disagree about the position of this limit, in part recognized in the route of the *Via Annia* [13], according to others, it should be searched further north [14]. Of the latter opinion were the scholars who, basing on the photo-interpretation of some aerial images, recognized in some anomalies visible in Ceggia and in the neighboring territories the result of hydraulic-agricultural interventions to be related to the nearby settlement of *Civitas Nova*/*Heraclia* (a center developed starting from the 7th century AD in a Byzantine lagoon environment) (Figure 2a), though not excluding a possible Renaissance origin [14,15]. The question, therefore, is now open. At the current state of knowledge, it is not possible to establish which sphere of influence falls this area or define its ancient environmental configuration.

Moving from previous considerations, we started to analyze some undefined artificial and natural evidence visible in multitemporal aerial images. Similar traces are identified in the Venetian plain [1,14], and their analysis could reveal important archeological infrastructures and sites [5,6,16,17]. More in general, these systems of anthropogenic traces are generally found in alluvial and coastal plains and often related to ancient settlements [18,19]. For this study, with the specific aim of defining the nature, chronology, and function of these undefined buried structures, in particular, we integrated multi-temporal aerial photo interpretation, geophysical measurements, core sample analysis, 14C dating, analysis of historical cartography, and archaeological survey [20].

From a geophysical point of view, the specific characteristics of the current area and the open questions related to wet soils suggested the application of FDEM and ERT methods to identify both natural and anthropogenic features visible from aerial and satellite images.

Electromagnetic methods in the frequency domain, known as FDEM or electromagnetic induction methods (EMI), are widely used for soil mapping in order to obtain a quick overview of the possible heterogeneity of a given conductive system [21–27]. One of the main variables that characterize soils, determining both vertical and horizontal variability within them, beyond the grain size and composition, is undoubtedly the water content. The measurement of the electrical conductivity of soils is, in fact, closely related and dependent on their water content/moisture content. The variation in the composition of a given soil, on the other hand, determines the ability to retain the moisture content di fferently. In general, apart from clayey soils, whose electrical response of high conductivity is inherent like the material, the di fferent grain size of the soil and the water content (saturated or unsaturated), determines the di fferent electrical response of the investigated system [28–32]. In particular, soils with smaller grain size will be more conductive than those with larger ones. The possibility of detecting fast, non-direct information on these characteristics with EM methods is, therefore, at the basis of their growing popularity in contexts where these data can contribute, for example, to the optimization of cultivation practices [21,33,34].

The main advantage o ffered by techniques able to indirectly measure the electrical conductivity of the investigated systems, such as low-frequency EM measurements (FDEM or EMI) and electrical resistivity tomography (ERT), is inherent in the ability to establish the spatial distribution and output of bodies characterized by a di fferent electrical conductivity (or its inverse resistivity) in an indirect way and with a di fferent degree of detail (resolution) and at di fferent depths. In particular, the strength of the FDEM technique lies undoubtedly in the speed of acquisition of the data of apparent conductivity of soils, made possible by the fact that direct contact with the measurement system is not necessary and therefore, the operator or the motorized vehicle carrying the instrument, connected to a GPS, can map quite quickly even large extensions giving general information on the heterogeneity of the system. The ERT method instead allows us to define in detail the exact depth and the relationship between the di fferent buried structures in the same system. The possibility to combine these two geophysical methods through comparable measurements allow to identify targets and systems of interest. This data integration is being increasingly used in the context of geomorphological and archaeological studies. This combination, in fact, allows to identify and describe both in a timely manner with FDEM methods natural structures related to the presence of water (e.g., paleo-environment), where the ERT method helps to better detail those and to recognize the presence of anthropogenic remains, whose exact identification requires a higher resolution than that o ffered by EM methods [35–43].

The more general objective of this study is to set up and sugges<sup>t</sup> a workflow reproducible in other contexts to allow the proper identification, valorization, and protection of complex archaeological landscapes similar to this here analyzed.

#### **2. Materials and Methods**

#### *2.1. Analysis of Historical Documents and Archaeological Field Data*

A review of the published and unpublished material kept at the Archive of the Soprintendenza Archeologia, Belle Arti e Paesaggio (SABAP) for the metropolitan area of Venice and the provinces of Padua, Belluno, and Treviso was firstly carried out.

On the other hand, between 2014 and 2017, new data was collected in the area [20,44] through systematic archaeological surface surveys [45–47].

This data was then gathered within a GIS platform, correctly structured to manage archaeological information; and the geomorphological picture, historical cartography, and aerial image traces were analyzed.

#### *2.2. Multitemporal Analysis of Historical Aerial Photographs*

Historical aerial photographs often represent a significant resource to identify the transformations that a ffected the landscapes during the times, changing their original or ancient asset. The visibility of interesting potential traces on the surface linked to the past by these supports is undoubtedly conditioned by many factors as, for example, the light, the vegetation, the ground humidity, the altitude of the flight, the land use, etc. [48,49].

For our purposes, we use these materials to:


The first step of the survey consisted of the identification of the aerial photographs preserved in the archives of various entities present on the national territory (Military Geographic Institute-IGM; Centre for Cartography of the Veneto Region-ReVen).

For the above reasons and purposes, the choice of the photos (Table 1) consider the number and quality of visible anomalies, and the temporal coverage, to allow a multitemporal analysis of the landscape. The most significant images, in this sense, result in the ReVen 1983, ReVen 1990, and GAI 1954 flights (Figures 3 and 4a), while the ReVen 1999 flight returned fewer details while confirming the presence of the most relevant evidence. Even the oldest images of the IGM 1937 flight reveal the anomalies but cover only a portion of the study area.


**Table 1.** List of aerial photographs analyzed for the study.

**Figure 3.** Area SE of Ceggia: zoom on some evident soil marks visible in historical zenithal photos of different years and a recent oblique photo. In yellow, the area subject to geophysical investigations, core drilling, and cleaning of the exposed sections. (**a**) Photo GAI 1954. (**b**) Photo ReVen 1983. (**c**) Photo ReVen 1990. (**d**) Oblique photo [5].

**Figure 4.** SE area of Ceggia. (**a**) System of traces visible in ReVen 1983 aerial frames. (**b**) Photorestitution on CTR 1:5000 of the signs detected through multitemporal analysis of historical aerial images (IGM 1937, GAI 1954, ReVen 1983 and Reven 1990) and new archaeological data deduced through surveys (2014–2017).

When digitalized, the frames of different flights have been geo-referenced in GIS in order to proceed with the identification of the anomalies and to compare with other data, (e.g., vector and raster types) available for the study area (geomorphological, archaeological, geophysical, cartographic, etc.).

The visible traces in the different temporal images were then divided into different layers, from the most recent (ReVen 1999) to the most dated (IGM 1937).

The information about the location, the source used to identify the anomaly, type of anomaly, and, when it possible, the interpretation and chronology have been associated.

#### *2.3. Multitemporal Analysis of Historical Cartography*

Historical cartography provides information on the shapes and conditions of the areas, their nature, but also on works planned realized or not realized. In our case, these instruments help to obtain information useful to interpret the traces visible in aerial images. These materials were mainly searched at the Archives of the Consorzio di Bonifica Veneto Orientale (CBVO) and the Archivio di Stato of Venice (ASVE), considering that the area of interest insists in a reclamation zone and that, since the Republic of Venice, it has been a ffected by operations aimed at regulating the waters. In the first case, project drawings, cadastral maps, and plans with sporadic indications about the use of defined sectors of the study area and on specific interventions concerning the hydrographic network were found [50]. The Archives of Venice, on the other hand, provide a part of the historical cartography make from the 16th century, by specialized institutes set up within the Republic of Venice to deal with the regulation of water, but also with reclamation, irrigation, and land culture. The same archive also contains maps of the Austrian Land Registry containing the description and use of the di fferent areas. A detailed representation of the nineteenth-century natural characteristics and anthropogenic structure of the area is provided by the map of the Second Military Survey and the Kriegskarte. Other significant detailed and geometrically correct restitution of the territory between 1890–1952 has been recovered by the cartography of the IGM, that has been mapping the entire Italian territory since the second half of the 19th century.

Di fferent cartographic products (Table 2) have then been inserted and georeferenced in a GIS project for the analysis and comparison with aerial photographs (Figure 5).


**Table 2.** List of historical cartography analyzed for the study.


**Table 2.** *Cont.*

**Figure 5.** Historical maps of the study area (in orange), and details of the route of the *Via Annia* (in yellow) and the Roman bridge (in black). (**a**) Drawing (1568) showing the sector to be reclaimed, i.e., the "retrazer palv", south of the Piavon canal (today's Canalat)-Archivio di Stato di Venezia, MP 42-55 (16); (**b**) Map (1639) depicting the marshes (in light yellow), ploughed and inhabited land (in pink), watercourses (in blue)-Archivio di Stato di Venezia, SEA, Piave, r.106, dis.15; (**c**) Detail of the map (1840) of the Santini Cadastre (in blue), superimposed on the CTR 1:10,000 (in black), which shows the relationship of the Canalat with the *Annia* Roman bridge; (**d**) Detail of a 1914 map showing the rectification of the Canalat bend in relation to the Roman bridge-Consorzio di Bonifica Veneto Orientale. The route of the *Via Annia*, the Canalat course and the location of the Roman bridge allow to identify the study area in the figures.

We used the geographical coordinates of the IGM tablets to identify clear ground control points (GCP) useful to geo-refer the historical cartography.

At the end of this process, we produced a multi-temporal sequence of images with a window of the transformations in the area of study from the mid-16th century to the present.

## *2.4. Geophysical Measurements*

In February 2015, we planned an FDEM field campaign to easily and quickly identify the main buried structures visible in the aerial and satellite images. Starting from these pieces of evidence, we selected a representative area, with some natural and artificial features, useful to collect field geophysical data. For the acquisition of FDEM data, we used a CMD electromagnetic conductivity meter (from GF Instruments s.r.o. Brno, Czech Republic) with a CMD-1 single-depth probe, setting high depth configuration, corresponding to −1.5 m full depth range. The resolution of the used instrument is 0.1 mS/m, where the accuracy corresponds to ±4% at 50 mS/m. The CMD-1 probe was manually moved on the ground in GPS continuous acquisition mode, connecting the instrument with a Trimble 5800 GPS receiver to register the right position of the acquisition lines. The FDEM data was collected on the field f.9 (Figure 6b) every meter, in a rectangular area, along parallel lines NE-SW oriented, for a total length of about 210 m. The entire width of the investigation area in the NW-SE direction was about 100 m. Due to the presence of small drainage channels, NE-SW oriented, the survey area has been divided into three main parts with a width of about 33 m for each (Figure 6). The conductivity values, during the processing, were converted in resistivity and then together plotted using the Surfer 10 Golden Software to obtain a map of the mean distribution of this parameter in the entire investigated area.

**Figure 6.** Geophysical measurements: (**a**) Image of an ERT line acquisition; (**b**) Localization of geophysical measurements: blue line (1) ERT, yellow area (2) FDEM; (**c**) scheme of ERT acquisition.

In the middle of the same field (f.9 in Figure 6), an electrical resistivity tomography (ERT) was acquired in the SW-NE direction. The electrical resistivity tomography (ERT) was performed by means of 15 lines, using for each line 48 electrodes spaced 0.5 m apart, overlapping line by line of 24 electrodes (roll-along) over a total length of 191.5 m. Direct and reciprocal measurements for each ERT line were acquired, by swapping current with potential electrodes, to estimate the errors in the dataset [51]. A di fference between di fferent cycle results (quality factor "Q") equal to 5% was imposed. The ERT inversion was performed using a regularized weighted least squares approach, according to the Occam's rule [52] as proposed by LaBrecque et al. [53]. The smoothness of the resistivity distribution here obtained strictly depends on the errors in each dataset.

According to Binley et al. [54], the best evaluation of errors might be obtained thanks to direct and reciprocal measurements, where these measurements shall be equal, providing the same resistance value. The possible deviation may be interpreted as an error estimate useful for the inversion. In this case study, the errors were smaller than 5%. The visualization of the inverted ERT lines was done using Surfer 10 Golden Software.

#### *2.5. Core Sample Analysis and Radiocarbon Dating*

Di fferent core samples were carried out using an Edelman combination-type hand auger. Sediment description included grain size, sedimentary structures, Munsell Charts color, presence, size and abundance of pedologenic nodules and mottles, paleontological content (shells, plant remains, wood, charcoal), archaeological content (sherds, bricks etc.). The data were compiled in stratigraphic logs and sections (Figures 7 and 8). An AMS radiocarbon date was obtained from a wood fragment (Table 3) that was collected from the undisturbed inner part of core CEG 7b and stored in aluminum foil after extraction in order to avoid contamination. The radiocarbon date was calibrated using software OxCal 4.3 [55] and applying the calibration curve IntCal13 [56].

**Figure 7.** Drilling. (**a**) Location ofthe cores: (1) andERT alignment (2). (**b**) Stratigraphic logs of core drillings performed in correspondence of paleochannels (4, 9, 8), ditches (7b), and main axes of the presumed stractures of ancient origin (9, 2a). Log symbols: (1) fragments of plant remains; (2) continental mollusks; (3) carbonatic nodules; (4) nodules of iron-manganese oxides; (5) charcoal; (6) bricks fragments. On the right, from top to bottom, photos of the works: (**c**) execution of a core; (**d**) sampling of a core; (**e**) documentation of stratigraphic succession.

**Figure 8.** Section across the study site.

**Table 3.** Details of the 14C dating of a wood sample from the core 7b (see Figure 7).


#### *2.6. Visual Evaluation of Exposed Section*

The exposed section of the southern side of the drainage channel that limits the field f.9, where geophysical prospecting and the cores drilling were realized, allowed us to analyze the correspondent section of the traces visible in the aerial photo (Figure 9a). Thanks to this activity, we documented the stratigraphic sequence [57] of the deposits and obtained detailed information about the different soil marks revealed in the aerial photos (e.g., potential road axes, drainage ditches, and paleochannel). The position of the sections of interest were identified thanks to the georeferenced images of the area stored in the GIS platform, where their coordinates were collected using a differential GPS and positioned on the field with several material references (i.e., pickets). The section is documented by photos, drawings, and stratigraphic units, interested up to 0.9 m depth from the ground level.

**Figure 9.** Cleaning of the exposed section. (**a**) Working area: (1) cores; (2) ERT line; (3) exposed sections. (**b**) Image documenting the field activity. (**c**) Stratigraphic sections detected in correspondence of some soil marks: the numbers correspond to the Stratigraphic Units (US) identified; the image shows the included angular US104 (hypothetical road/street base).
