Shallow Offshore Geophysical Prospection of Archaeological Sites in Eastern Mediterranean

Geophysical prospecting methods have been extensively used to outline buried antiquities in terrestrial sites. Despite the frequent application of these mapping and imaging approaches for the detection of archaeological relics in deep-water marine environments (e.g., shipwrecks), the aforementioned processes have minimal contribution when it comes to understanding the dynamics of the past in coastal and shallow aquatic archaeological sites. This work explores the possibilities of multicomponent geophysical techniques in revealing antiquities that have been submerged in diverse shallow coastal marine environments in the eastern Mediterranean. A group of four sites in Greece (Agioi Theodoroi, Olous, Lambayanna) and Cyprus (Pafos) spanning from prehistory to Roman times were chosen as test sites to validate the efficiency of electrical resistivity tomography, magnetic gradiometry, and ground penetrating radar methods. The comprehensive analysis of the geophysical data completed the picture for the hidden archeological elements in all the sites. The results manifest the significance and the potential of these methods for documenting and understanding the complex archaeological sites encountered in the Mediterranean. In view of climate change and the risks related to future sea level rise and erosion of low-level coastal areas, the results of this work could be integrated in a strategic framework to develop an effective interdisciplinary research model that can be applied to similar shallow water archaeological surveys, thus substantially contributing towards cultural resources management.


Introduction
The study and documentation of Cultural Heritage (CH) is rapidly developing, at times favoring or even pushing the progress of technologies and equipment in place for the digitalization of items and landscapes. Nowadays, it is possible to document in high detail our surroundings, making a step further to "look" beneath the soil and produce images of buried structures to assist and complement archaeological research in a timely and costly efficient manner [1].
Underwater investigation, although still a younger branch of archaeology, was the leader of rapid and important development as well. The technological and methodological innovations allowed it to pass from an "amateur" activity to a proper scientific field of research with substantial contribution in the domains of past commerce and coast exploitations. In this sense, most of the efforts so far have been placed in deep water context explorations, where specific factors that do not apply in other archaeological research-like pressure, visibility, oxygen availability-need to be taken into account for the successful implementation of the research activity [2].
However, the coastal and shallow underwater area is still considered as a "gray zone", where different factors create a unique, constantly evolving and challenging environment, sharing elements from both deep-water investigations and dry land surveys. Indeed, the littoral area with its associated archaeology is often easily accessible by experts and the

Materials and Methods
The research was undertaken in hierarchical and subsequent modules, where practical applications were mixed with testing, improvement, and optimization. The first phase consisted of proper documentation and digitization of legacy data with the objective to gather the existing historical and scientific data from previous archaeological references, fieldwork projects, and related publications. The digitization and processing of topographical plans, site maps, and the archaeological data from published material into Geographical Information Systems (GIS) using ArcGIS version 10.4 aimed at facilitating the integrated analysis of every different archaeological site. At the same time, a Real-Time Kinematic Global Navigation Satellite System (RTK-GNSS) employing two Javad Triumph-1 sensors working as base and rover systems were used to map and document the in situ visible onshore and submerged archaeological structures (walls, pavements, buildings, corridors etc.) with sub-centimeter accuracy and up to a manageable depth of the littoral zone. This module resulted in the creation of fully geo-referenced site plans based on fixed ground control points, including all archaeological features of the different sites.
In parallel to this, Remotely Piloted Aerial Systems (RPAS) using a kite and a Phantom-3 drone were used to extensively document the sites by means of digital photogrammetry and accurate image georeferencing. The compiled digital elevation model for each site was based on the digital photogrammetry and was instrumental for the documentation of their status, but also as a calibration bench for the subsequent multicomponent geophysical mapping. The Agisoft Metashape software was used for the photogrammetric data processing.
The geophysical survey was concentrated on the coastal and the shallow marine part of the sites in an effort to complete the picture of the onshore and offshore hidden built environment. The strategy of the geophysical survey was appropriately adjusted to meet the challenges and questions that need to be answered in each of the sites. Thus, a multicomponent geophysical imaging approach integrating different methods was efficiently organized. More specifically, ground penetrating radar was employed on the coast to map the horizontal extension and dimension of the architectural relics up to the depth of 1.5-

Materials and Methods
The research was undertaken in hierarchical and subsequent modules, where practical applications were mixed with testing, improvement, and optimization. The first phase consisted of proper documentation and digitization of legacy data with the objective to gather the existing historical and scientific data from previous archaeological references, fieldwork projects, and related publications. The digitization and processing of topographical plans, site maps, and the archaeological data from published material into Geographical Information Systems (GIS) using ArcGIS version 10.4 aimed at facilitating the integrated analysis of every different archaeological site. At the same time, a Real-Time Kinematic Global Navigation Satellite System (RTK-GNSS) employing two Javad Triumph-1 sensors working as base and rover systems were used to map and document the in situ visible onshore and submerged archaeological structures (walls, pavements, buildings, corridors etc.) with sub-centimeter accuracy and up to a manageable depth of the littoral zone. This module resulted in the creation of fully geo-referenced site plans based on fixed ground control points, including all archaeological features of the different sites.
In parallel to this, Remotely Piloted Aerial Systems (RPAS) using a kite and a Phantom-3 drone were used to extensively document the sites by means of digital photogrammetry and accurate image georeferencing. The compiled digital elevation model for each site was based on the digital photogrammetry and was instrumental for the documentation of their status, but also as a calibration bench for the subsequent multicomponent geophysical mapping. The Agisoft Metashape software was used for the photogrammetric data processing.
The geophysical survey was concentrated on the coastal and the shallow marine part of the sites in an effort to complete the picture of the onshore and offshore hidden built environment. The strategy of the geophysical survey was appropriately adjusted to meet the challenges and questions that need to be answered in each of the sites. Thus, a multicomponent geophysical imaging approach integrating different methods was efficiently organized. More specifically, ground penetrating radar was employed on the coast to map the horizontal extension and dimension of the architectural relics up to the depth of 1.5-2.0 m below the ground surface. Magnetic gradiometry covered extensive offshore sections of the sites in an effort to compile interpretation maps showing the distribution of architectural features located beneath the sea bottom up to the depth of about one meter. Dynamic and/or static submerged 3-D Electrical Resistivity Tomography (ERT) provided stratigraphic sections up to the depth of~3-4 m below the seabed, including the vertical dimension in the reconstruction of the submerged buried cultural material. The conductivity of the seawater was accurately measured with a high precision/resolution handheld conductivity meter. The bathymetry model for the very shallow part was compiled by moving the rover GNSS unit along parallel transects in the sea and placing the pole in different points on the seabed.
All the spatial data (digitized maps, plans, aerial images, bathymetric models, geophysical maps, and interpretation) were imported in a GIS platform for further analysis and processing. This enabled the compilation of new thematic maps, the identification of patterns in the built environment, the extraction of spatial information (e.g., dimensions, alignment) of the structures. Collected and analyzed data were condensed in a comprehensive archaeological interpretation with the aim to reconstruct possible chronological phases, potential continuity of use, and space organization of identified buried and/or submerged structures.

Agioi Theodoroi, Crete, Greece
The coastal archaeological area of Agioi Theodoroi is located on the northern coast the island of Crete (Figure 1). A standard photo camera was mounted on a kite to create a high-resolution aerial image of the site and map the current status of the visible relics on land and those submerged in the sea. The photographs were taken at an altitude of between 50 and 100 m, estimated through the photogrammetric processing. More than 170 photos were selected for inclusion in the final photogrammetric model and combined into a composite orthophoto using commercial software. The orthophoto was georectified using ground control points mapped with a differential GNSS unit with a horizontal accuracy of less than 50 cm (Figure 2a). The visible submerged archaeological walls (Figure 2e) were mapped with the GNSS system and superimposed on the aerial photo as red polylines (Figure 2a).
The archaeologist Spiros Marinatos completed the first systematic survey and excavation of the site in the early 20th century [14] and brought to light a late Minoan prehistoric seaside building as well as wall constructions that continue towards the sea. The older archaeological plan of the site, including the excavated and visible architectural relics from the 1924 excavation campaign, was rectified using common points with the aerial image from the quarry to the north and the visible remains of the small building on the sand beach ( Figure 2b). The rectification of this plan, employing an affine transformation, included positioning errors of~5% due to the limited common points between the aerial image and the respective plan. However, these efforts comprised the most efficient approach to document the excavated building complexes along with visible structures on the modern landscape ( Figure 2a). A rocky outcrop to the north with three consecutive islets protects the site from the respective winds ( Figure 2c). The specific islets show strong traces and indications of quarrying activities. To the west of this outcrop (Figure 2d), there is a rectangular carving that is divided into two parts by a strip of uncut rock. The excavator considered this as the best natural jetty in the area and the large carving as a slipway or shipshed. Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 18 The Pole-Dipole electrode array in forward and reverse mode and an electrode spacing of one meter were used to collect the apparent resistivity tomographic data along all the profiles. For the marine tomographic sections, the multinode cable was submerged and laid on the bottom of the sea during data collection. The cable was fixed with two anchors at both ends during the measurement to avoid unnecessary movement due to the waves. The seawater conductivity was measured to 0.2 Ohm-m with a high sensitivity and accuracy portable conductivity meter. The survey pattern of colleting a number of parallel profiles with line spacing three times larger than the basic electrode spacing in the marine grid prohibited the use of a 3-D inversion approach to reconstruct the resistivity structure below the sea-bottom. Thus, these profiles were processed only with a 2-D robust inversion algorithm [15]. The resistivity of the sea water column and the depth information for each electrode location from the sea water surface were also incorporated into The Pole-Dipole electrode array in forward and reverse mode and an electrode spacing of one meter were used to collect the apparent resistivity tomographic data along all the profiles. For the marine tomographic sections, the multinode cable was submerged and laid on the bottom of the sea during data collection. The cable was fixed with two anchors at both ends during the measurement to avoid unnecessary movement due to the waves. The seawater conductivity was measured to 0.2 Ohm-m with a high sensitivity and accuracy portable conductivity meter. The survey pattern of colleting a number of parallel profiles with line spacing three times larger than the basic electrode spacing in the marine grid prohibited the use of a 3-D inversion approach to reconstruct the resistivity structure below the sea-bottom. Thus, these profiles were processed only with a 2-D robust inversion algorithm [15]. The resistivity of the sea water column and the depth information for each electrode location from the sea water surface were also incorporated into the inversion model to account for the effect of the highly conductive saline water layer. On the other hand, the land profiles were combined and processed with a 3-D inversion algorithm [16,17] to produce horizontal resistivity slices in different depths.
The 2-D inverted resistivity models along Lines D01 to D14 show a two-layer stratigraphy up to the depth of about 8 m below the free surface of the sea. Below the water layer, that has an average thickness of 0.8 m, there is a relative conductive layer with resistivity less than 1 Ohm-m and average thickness of about 2 m, which can be correlated with the Pleistocene to Holocene coastal sands. A more resistive layer (1.5-7 Ohm-m) represents the geological background of the area composed of the upper Miocene marly limestone. Within the layer of the coastal sands, the resistivity sections outline isolated resistive anomalies (outlined by the white circles in Figure 3), which seem like outcrops of the bottom layer reaching up to the bottom of the sea. These resistive targets, like those for example in Line D03, correlate well with visible submerged walls. This in situ acted as a guide for isolating resistive targets with similar signatures that were attributed to "unseen" ancient walls, thus completing in this way the picture regarding unknown submerged archaeological structures ( Figure 3).
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 18 the inversion model to account for the effect of the highly conductive saline water layer. On the other hand, the land profiles were combined and processed with a 3-D inversion algorithm [16,17] to produce horizontal resistivity slices in different depths. The 2-D inverted resistivity models along Lines D01 to D14 show a two-layer stratigraphy up to the depth of about 8 m below the free surface of the sea. Below the water layer, that has an average thickness of 0.8 m, there is a relative conductive layer with resistivity less than 1 Ohm-m and average thickness of about 2 m, which can be correlated with the Pleistocene to Holocene coastal sands. A more resistive layer (1.5-7 Ohm-m) represents the geological background of the area composed of the upper Miocene marly limestone. Within the layer of the coastal sands, the resistivity sections outline isolated resistive anomalies (outlined by the white circles in Figure 3), which seem like outcrops of the bottom layer reaching up to the bottom of the sea. These resistive targets, like those for example in Line D03, correlate well with visible submerged walls. This in situ acted as a guide for isolating resistive targets with similar signatures that were attributed to "unseen" ancient walls, thus completing in this way the picture regarding unknown submerged archaeological structures ( Figure 3). The small rectangular grid on the coast shows a high inhomogeneous and conductive background material with an average resistivity of about 3 Ohm-m. The resistivity values lower than the background signify the sandy superficial layer that is saturated with the saline seawater. High resistivity values appear at the southern part of the grid forming a compact region that reaches the depth of 2 m below the surface. This resistive region is confined within the boundaries of the visible walls of the archaeological building. Although resistivity depth slices managed to outline the general shape of the building, the ERT survey in this grid failed to delineate any other structural details due to the survey parameters used to collect the tomographic data using a line separation twice the interelectrode distance along the profiles (Figure 4). The small rectangular grid on the coast shows a high inhomogeneous and conductive background material with an average resistivity of about 3 Ohm-m. The resistivity values lower than the background signify the sandy superficial layer that is saturated with the saline seawater. High resistivity values appear at the southern part of the grid forming a compact region that reaches the depth of 2 m below the surface. This resistive region is confined within the boundaries of the visible walls of the archaeological building. Although resistivity depth slices managed to outline the general shape of the building, the ERT survey in this grid failed to delineate any other structural details due to the survey parameters used to collect the tomographic data using a line separation twice the inter-electrode distance along the profiles (Figure 4). Line B1 shows high resistivity values on the 10th, 17th, 26th, 29th, 31st, 38th, and 40th meter from the beginning of the profile to the east. The black rectangular areas in Figure  5 show traces of small walls located on the seabed and seem to extend to a depth of about 0.7-1.0 m below the sea bottom. Line B2 outlines similar small resistive targets in the superficial layer below the sea bottom showing the horizontal extent of these walls between the two lines ( Figure 5). The inclined and relativity more compact resistivity surface characterizes the western part of both profiles from the 27th meter until the end of the lines. This surface is presented as a black dashed dipping line in Figure 5 and shows the rocky bedrock that is dipping to the east with a slope of about 26 degrees. The morphology, the resistive signature, and the slope of the rocky bedrocks are concrete indications to support the hypothesis for the function of this part as a slipway or shipshed ( Figure 5). Line B1 shows high resistivity values on the 10th, 17th, 26th, 29th, 31st, 38th, and 40th meter from the beginning of the profile to the east. The black rectangular areas in Figure 5 show traces of small walls located on the seabed and seem to extend to a depth of about 0.7-1.0 m below the sea bottom. Line B2 outlines similar small resistive targets in the superficial layer below the sea bottom showing the horizontal extent of these walls between the two lines ( Figure 5). The inclined and relativity more compact resistivity surface characterizes the western part of both profiles from the 27th meter until the end of the lines. This surface is presented as a black dashed dipping line in Figure 5 and shows the rocky bedrock that is dipping to the east with a slope of about 26 degrees. The morphology, the resistive signature, and the slope of the rocky bedrocks are concrete indications to support the hypothesis for the function of this part as a slipway or shipshed ( Figure 5).

Ancient Olous, Crete, Greece
Ancient Olous was an important center and a strategic port located on the northern coast of Crete. According to the archaeological evidence, the area exhibited a continuous habitation spanning from the Minoan to the Venetian period. The most relevant architectural constructions, though, refer to the Hellenistic, Roman, and Early Byzantine periods. The city was a flourishing harbor town with a sanctuary, a necropolis, and its own coinage since its earlier phases. Moreover, the epigraphic evidence from the area indicates that the Hellenistic Olous had been claimed both from the Ptolemies and the independent state of Rhodes, probably because it represented a strategic place on the Aegean Sea routes [18,19]. Because of both the eustatic sea level change and the tectonic depression of the area, most of the remains of the ancient structures that are attributable to the main archaeological occupation of Olous are now either submerged or covered with sand of the seabed.
A different ERT survey strategy was followed to map the rectangular section of 41 m by 20 m on the very shallow part where the water depth was less than one meter, in order to validate the efficiency of the floating survey mode in actual conditions. A dense network of parallel ERT lines were completed along a single direction and one-meter separation. The custom made multinode cable had 13 sensors with an interval of one meter and was kept on the surface of the sea with a number of floaters. The 10-channel electrical resistivity instrument, an external battery, and a toughpad were appropriately connected and placed inside a floating apparatus for protection. The integrated system was floated slowly along the parallel transects to complete the tomographic measurements. A dipoledipole array with one transmitting dipole close to the apparatus and 10 receiving dipoles with gradually increasing distance from the current electrode pair was used to capture and store the underwater apparent resistivity readings through commercial marine management software. The resistivity of the sea water column (0.17 Ohm-m) constrained the 3-D processing (Figure 6a).

Ancient Olous, Crete, Greece
Ancient Olous was an important center and a strategic port located on the northern coast of Crete. According to the archaeological evidence, the area exhibited a continuous habitation spanning from the Minoan to the Venetian period. The most relevant architectural constructions, though, refer to the Hellenistic, Roman, and Early Byzantine periods. The city was a flourishing harbor town with a sanctuary, a necropolis, and its own coinage since its earlier phases. Moreover, the epigraphic evidence from the area indicates that the Hellenistic Olous had been claimed both from the Ptolemies and the independent state of Rhodes, probably because it represented a strategic place on the Aegean Sea routes [18,19]. Because of both the eustatic sea level change and the tectonic depression of the area, most of the remains of the ancient structures that are attributable to the main archaeological occupation of Olous are now either submerged or covered with sand of the seabed.
A different ERT survey strategy was followed to map the rectangular section of 41 m by 20 m on the very shallow part where the water depth was less than one meter, in order to validate the efficiency of the floating survey mode in actual conditions. A dense network of parallel ERT lines were completed along a single direction and one-meter separation. The custom made multinode cable had 13 sensors with an interval of one meter and was kept on the surface of the sea with a number of floaters. The 10-channel electrical resistivity instrument, an external battery, and a toughpad were appropriately connected and placed inside a floating apparatus for protection. The integrated system was floated slowly along the parallel transects to complete the tomographic measurements. A dipole-dipole array with one transmitting dipole close to the apparatus and 10 receiving dipoles with gradually increasing distance from the current electrode pair was used to capture and store the underwater apparent resistivity readings through commercial marine management software. The resistivity of the sea water column (0.17 Ohm-m) constrained the 3-D processing (Figure 6a).
The horizontal depth slices extracted from the 3-D resistivity inversion model clearly outline the continuation of the wall seen in the aerial image running for a distance of 20 m towards the north and up the coastline. The resistivity imaging also outlines at least two parallel east-west trending walls, that are oriented vertically to the previous south-north structure, reaching the depth of 0.6 m below the seabed. The compact resistive region at the east, which is clearly observed up to the depth of 1.2 m, is probably related to a collapsed structure with a relatively poor preservation level (Figure 6b). The horizontal depth slices extracted from the 3-D resistivity inversion model clearly outline the continuation of the wall seen in the aerial image running for a distance of 20 m towards the north and up the coastline. The resistivity imaging also outlines at least two parallel east-west trending walls, that are oriented vertically to the previous southnorth structure, reaching the depth of 0.6 m below the seabed. The compact resistive region at the east, which is clearly observed up to the depth of 1.2 m, is probably related to a collapsed structure with a relatively poor preservation level (Figure 6b).
The results from the experimental ERT survey in Olous during 2016 comprised the first systematic effort in mapping non-visible submerged archaeological relics in shallow The results from the experimental ERT survey in Olous during 2016 comprised the first systematic effort in mapping non-visible submerged archaeological relics in shallow marine environments where the water thickness is less than one meter. Based on the previous experience with extensive ERT numerical modeling [20], the method was successfully transferred in real conditions and successfully managed to reveal the hidden submerged ancient walls as resistive targets within the highly conductive marine context.

Lambayanna, Peloponnese, Greece
The submerged settlement in the bay of Lambayanna (Figure 7a) is dated to the first half of the 3rd millennium BCE (Early Bronze Age II). The initial underwater archaeological survey in 2015 indicated the expansion of the settlement within a marine area was of covering an area of at least 12,000 m 2 . About 6000 artifacts consisting of ceramics, stone tools, and obsidian blades were recovered from the sea bottom. The in situ survey and preliminary underwater archaeological mapping revealed scattered walls and stones all over the bay in a water depth ranging from 1 m to 3 m. At least three massive U-shaped aligned structures and stretches of fortification walls were found [21]. The underwater excavation activities undertaken within two stratigraphic soundings brought to light ceramics dated to the Early Bronze Age, meaning the late 4th or early 3rd millennium BCE [22]. The central part is occupied by a relatively more disturbed magnetic signal in relation to the southern section, which is related to natural processes and the material that is deposited in the sea from the local seasonal stream that flows from the east along the dirty road leading to the beach. However, there are some evidences of cultural material represented by both negative (M7, M9, M10) and positive (M8, M11) anomalies that share similar SW-NE orientation. The structural complex M8 occupies an area of 11 × 16 m and is divided into at least three compartments. The S-N linear segments M9 and M10 outline potential walls that probably belong to the same structural phase due to their similar orientation and magnetic signature. A couple of positive anomalies (M11) that slightly deviate from the wider orientation are also outlined ( Figure 8). The geophysical survey in the bay was concentrated in the area outlined by the red polygon outlines in Figure 7a. A series of eight magnetometers 0.25 m apart were mounted on a non-metallic frame and connected to a GNSS rover station though an associated electronic box (Figure 7b). The whole system was mounted on a custom-made floating plastic apparatus, which also facilitated the vertical adjustment of the frame in order to adjust the position of the sensors related to the seabed. This adjustment was very important in cases of submerged obstacles (e.g., large stones) that could damage the sensors in cases of sudden impact. Specific managing software controlled the data acquisition along almost south-north tracks covering more than 9000 m 2 of the shallow marine section. The magnetic data were processed following a systematic flowchart including the removal of extreme magnetic values and overlapping points. The median value was subtracted from the magnetic values along each sensor transect in order to remove the stripes, due to spatially close or overlapping transects [23]. A cubic spline interpolation was used to compute new grid nodes every 0.1 m in both directions and, finally, a Gaussian low pass filter was employed to remove the high frequency noise, thus resulting in a smoother grid.
The central part is occupied by a relatively more disturbed magnetic signal in relation to the southern section, which is related to natural processes and the material that is deposited in the sea from the local seasonal stream that flows from the east along the dirty road leading to the beach. However, there are some evidences of cultural material represented by both negative (M7, M9, M10) and positive (M8, M11) anomalies that share similar SW-NE orientation. The structural complex M8 occupies an area of 11 × 16 m and is divided into at least three compartments. The S-N linear segments M9 and M10 outline potential walls that probably belong to the same structural phase due to their similar orientation and magnetic signature. A couple of positive anomalies (M11) that slightly deviate from the wider orientation are also outlined (Figure 8). The northern part shows extremely high magnetic values related to metallic anchors located in the sea bottom, which are related to the use of this part as a small mole for local boats. The impact of the natural processes is also pronounced in this section of the site since the magnetic field has been severely disturbed due to the inhomogeneous material carried and deposited in the sea bottom by the adjacent seasonal stream. Within this extremely noisy environment, the gradiometer map seems to outline a poorly preserved rectangular structure (M13) with dimensions 19 × 11 m as negative values and a smaller rectilinear feature (M12) with a positive signature (Figure 8).
The Ground Penetrating Radar (GPR) covered an area of more than 1500 m 2 along The northern part shows extremely high magnetic values related to metallic anchors located in the sea bottom, which are related to the use of this part as a small mole for local boats. The impact of the natural processes is also pronounced in this section of the site since the magnetic field has been severely disturbed due to the inhomogeneous material carried and deposited in the sea bottom by the adjacent seasonal stream. Within this extremely noisy environment, the gradiometer map seems to outline a poorly preserved rectangular structure (M13) with dimensions 19 × 11 m as negative values and a smaller rectilinear feature (M12) with a positive signature (Figure 8).
The Ground Penetrating Radar (GPR) covered an area of more than 1500 m 2 along parallel transects every 0.5 m at the southern part of the bay using the 250 MHz antenna (Figure 7c). Standard processing routines were employed to enhance the GPR signal including trace reposition, time zero correction, dewow filter background subtraction, and exponential compensation gain. The velocity of the electromagnetic signal was calculated through the fitting of a hyperbola in respective reflections which appeared in the radargrams. Then, a Hilbert transform was applied to calculate the instantaneous amplitude and extract depth slices based on the calculated velocity of 0.1 m/ns.
Despite the saline nature and the increased conductivity of the background soil along the coast, the GPR slice from depth 0.2-0.3 m shows an area with bulk strong reflections (G1) which are either related to the concentration of boulders below the surface or to collapsed structural remains. The shape and linearity of the strong reflection G2 shows positive indication of a small rectangular building with dimensions 5 × 5 m. The two parallel linear segments (G3) further to the north that have comparable length (3.5 m) probably outline the location of small walls in depths less than 0.5 m. The group of linear strong reflections in G4 seem to belong to the same structural phase indicating isolated walls with similar length (~3-3.5 m) (Figure 9a). The central region shows a linear reflection (G5) at least 9 m long, probably related to architectural buried remains. Further to the NE of G5, the GPR slice from 0.5-0.6 m outlines a series of parallel linear segments most probably related to wall remains. The NE corner of the section is occupied by strong reflections that do not exceed the depth of 0.6 m, which could be attributed to cultural material (Figure 9b).
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 exponential compensation gain. The velocity of the electromagnetic signal was calculated through the fitting of a hyperbola in respective reflections which appeared in the radargrams. Then, a Hilbert transform was applied to calculate the instantaneous amplitude and extract depth slices based on the calculated velocity of 0.1 m/ns. Despite the saline nature and the increased conductivity of the background soil along the coast, the GPR slice from depth 0.2-0.3 m shows an area with bulk strong reflections (G1) which are either related to the concentration of boulders below the surface or to collapsed structural remains. The shape and linearity of the strong reflection G2 shows positive indication of a small rectangular building with dimensions 5 × 5 m. The two parallel linear segments (G3) further to the north that have comparable length (3.5 m) probably outline the location of small walls in depths less than 0.5 m. The group of linear strong reflections in G4 seem to belong to the same structural phase indicating isolated walls with similar length (~3-3.5 m) (Figure 9a). The central region shows a linear reflection (G5) at least 9 m long, probably related to architectural buried remains. Further to the NE of G5, the GPR slice from 0.5-0.6 m outlines a series of parallel linear segments most probably related to wall remains. The NE corner of the section is occupied by strong reflections that do not exceed the depth of 0.6 m, which could be attributed to cultural material (Figure 9b).

Pafos, Cyprus
The purpose of the geophysical survey in Pafos was to map the physical properties of the subsoil below the sea surface and to locate possible archaeological remains related to port facilities of the Hellenistic and Roman periods concerning stone masonry structures (Figure 10d). The geophysical survey covered a small rectangular grid with dimensions 15 m by 10 m using the 3-D resistivity tomography method (Figure 10a). The four

Pafos, Cyprus
The purpose of the geophysical survey in Pafos was to map the physical properties of the subsoil below the sea surface and to locate possible archaeological remains related to port facilities of the Hellenistic and Roman periods concerning stone masonry structures (Figure 10d). The geophysical survey covered a small rectangular grid with dimensions 15 m by 10 m using the 3-D resistivity tomography method (Figure 10a). The four corners of the survey area were stacked using RTK-GNSS, enabling high positioning accuracy. Multiple ERT parallel lines were developed in two vertical directions (x, y) in order to increase the density of the acquired information (Figure 10b). Parallel to the "x" direction 21 lines (0.5 m apart with 31 electrodes per line) and parallel to the "y" direction 16 lines (1.0 m apart with 21 electrodes per line) were installed. The static submerged mode was employed for the survey, where the cable was sunk and laid on the bottom of the seabed during data collection. The submerged cable was kept fixed with two anchors at both ends during the measurement to avoid unnecessary movement due to the waves. The Pole-Dipole and Gradient arrays were used for each line in both directions to collect the tomographic data. The specific survey pattern was selected based on previous research efforts signifying the advantage of combining submerged tomographic data that have been collected with multiple electrode arrays, in order to increase the resolving power of the final 3-D resistivity models below the seabed [20]. corners of the survey area were stacked using RTK-GNSS, enabling high positioning accuracy. Multiple ERT parallel lines were developed in two vertical directions (x, y) in order to increase the density of the acquired information (Figure 10b). Parallel to the "x" direction 21 lines (0.5 m apart with 31 electrodes per line) and parallel to the "y" direction 16 lines (1.0 m apart with 21 electrodes per line) were installed. The static submerged mode was employed for the survey, where the cable was sunk and laid on the bottom of the seabed during data collection. The submerged cable was kept fixed with two anchors at both ends during the measurement to avoid unnecessary movement due to the waves. The Pole-Dipole and Gradient arrays were used for each line in both directions to collect the tomographic data. The specific survey pattern was selected based on previous research efforts signifying the advantage of combining submerged tomographic data that have been collected with multiple electrode arrays, in order to increase the resolving power of the final 3-D resistivity models below the seabed [20].
The pole with the RTK-GNSS unit was placed on the position of every submerged probe of the cable to measure the respective depth. Then a kriging algorithm was used to compile the bathymetry model that shows a depth variation between 0.3 m to 2.3 m (Figure 10c). The conductivity value of the water was measured with handheld equipment in various positions and an average value of 0.19 Ohm-m was calculated. A 10-channel resistivity meter, an external battery, and a laptop were used to inject the current below the sea bottom through the sensors based on a predefined electrode protocol and measuring the resulting potential. The survey lines from both directions and electrode configurations were combined to a single file and a 3-D robust inversion algorithm [13] was used to reconstruct the resistivity variation below the sea bottom incorporating the resistivity and the thickness of the water layer in the inversion procedure. The combined interpretation of the resistivity depth slices outlined three different resistive areas of interest (R1, R2, and R3). The approximate depth of the targets is esti- The pole with the RTK-GNSS unit was placed on the position of every submerged probe of the cable to measure the respective depth. Then a kriging algorithm was used to compile the bathymetry model that shows a depth variation between 0.3 m to 2.3 m (Figure 10c). The conductivity value of the water was measured with handheld equipment in various positions and an average value of 0.19 Ohm-m was calculated. A 10-channel resistivity meter, an external battery, and a laptop were used to inject the current below the sea bottom through the sensors based on a predefined electrode protocol and measuring the resulting potential. The survey lines from both directions and electrode configurations were combined to a single file and a 3-D robust inversion algorithm [13] was used to reconstruct the resistivity variation below the sea bottom incorporating the resistivity and the thickness of the water layer in the inversion procedure.
The combined interpretation of the resistivity depth slices outlined three different resistive areas of interest (R1, R2, and R3). The approximate depth of the targets is estimated at 0.57-1 m below the seabed. Target R1a is parallel to the harbor wall with a direction SW to NE. The average width of the target is less than 0.5 m. Another smaller target attached to the aforementioned target (R1b) is located vertical to the main axis of R1a. Target R3 is located on the NE side of the area, vertical to the main axis of R1a with an average estimated width approximately 1 m. The rectangular structure R3b is attached to R3 and the SW side of the survey grid outlines where the R4 is located (Figure 11). mated at 0.57-1 m below the seabed. Target R1a is parallel to the harbor wall with a direction SW to NE. The average width of the target is less than 0.5 m. Another smaller target attached to the aforementioned target (R1b) is located vertical to the main axis of R1a. Target R3 is located on the NE side of the area, vertical to the main axis of R1a with an average estimated width approximately 1 m. The rectangular structure R3b is attached to R3 and the SW side of the survey grid outlines where the R4 is located (Figure 11).
The overall results in the port show strong evidence for the location of architectural structures related to buildings and structures that are submerged up to 2 m below the seabed. Within this layer, the geophysical survey was able to map potential building walls, which are in a relatively good conservation level, based on the strong and discrete resistive signature in the tomographic images. Thus, the results of this geophysical survey can form the basis for consolidating specific actions by the archaeological service to protect the archaeological material from future contraction works in the port.

Discussion and Conclusions
Low-level coastal zones are highly dynamic systems showing extensive complexity due to the effect of geological and climate change factors that act over different spatial and time scales [24]. The continuous action of wind and waves, the sediment movement, and associated sea level alterations cause a constant and progressive change on the coasts. The rise of the global sea level over the last centuries, in combination with the relevant projected rise of about 1 m by 2100 [25], will inevitably trigger severe flooding incidents and soil erosion along the shores [26,27].
These natural factors integrated with associated pressures from anthropogenic activities formulate a constantly evolving environment, rendering the littoral zones fragile and instable. Inevitably, the cultural resources lying in these zones face risks and threats for The overall results in the port show strong evidence for the location of architectural structures related to buildings and structures that are submerged up to 2 m below the seabed. Within this layer, the geophysical survey was able to map potential building walls, which are in a relatively good conservation level, based on the strong and discrete resistive signature in the tomographic images. Thus, the results of this geophysical survey can form the basis for consolidating specific actions by the archaeological service to protect the archaeological material from future contraction works in the port.

Discussion and Conclusions
Low-level coastal zones are highly dynamic systems showing extensive complexity due to the effect of geological and climate change factors that act over different spatial and time scales [24]. The continuous action of wind and waves, the sediment movement, and associated sea level alterations cause a constant and progressive change on the coasts. The rise of the global sea level over the last centuries, in combination with the relevant projected rise of about 1 m by 2100 [25], will inevitably trigger severe flooding incidents and soil erosion along the shores [26,27].
These natural factors integrated with associated pressures from anthropogenic activities formulate a constantly evolving environment, rendering the littoral zones fragile and instable. Inevitably, the cultural resources lying in these zones face risks and threats for deterioration and destruction. This generates the necessity to plan, organize, and implement actions to document and preserve the coastal and shallow submerged cultural material, assuring its accessibility and protection for the future generations [28].
The Minoan coastal area of Agioi Theodoroi was subject to a coastal and marine ERT survey. The rectification of the older topographic and excavation plan of the site on the aerial image, in combination with the geophysical results, verified the continuation of the visible onshore and submerged relics below the sea bottom, completing the picture of the structured environment and enhancing the hypothesis for the existence of a slipway or shipshed. The experimental floating dynamic ERT survey on a small grid in ancient Olous outlined a long linear resistive wall, which was crossed with at least three new vertical segments that were not visible with visual inspection. The coastal and marine geophysical campaign in Lambayanna provided promising results in the sections under investigation. Independently of modern or natural interventions to the site, the high-resolution images were able to indicate a number of candidate targets, thus completing to a certain degree the picture for the built environment that is buried below the seabed. The geophysical imaging survey used in the port of Pafos provided solid evidence on the existence of architectural structures related to building walls with a relatively good conservation level.
The successful results from different coastal and shallow submerged archaeological sites in the Mediterranean addressed the main objective of this work, proving the efficiency of the magnetic gradiometry, electrical resistivity tomography, and ground penetrating radar methods to reconstruct the built archaeological environment that lies a few meters below the seabed.
The specific framework generates the need to implement diverse geophysical imaging methods, integrated with complementary geoinformation techniques, to reconstruct the landscape and handle specific archaeological queries concerning the mapping of shallow offshore archaeological sites. Specifically, magnetic gradiometry was showed to be robust in outlining the horizontal extent of archaeological structural materials (e.g., walls, buildings) that lie in the sandy layer below the seabed. Two-dimensional electrical resistivity tomography can reconstruct the vertical stratigraphy below the sea bottom, resolving at the same time shallow submerged compact archaeological structures. Three-dimensional electrical resistivity modes, deployed either on the coast or in the shallow aquatic environment, were proven extremely efficient in delineating resistive cultural relics in a complete 3-D space. Even ground penetrating radar seems to extract useful archaeological outcomes, despite the conductive background soil matrix because of the increased salinity along the coastline.
The scientific originality of this work lies in extending the efficient employment of geophysical methods from terrestrial to ultra-shallow aquatic environments. This transition leads to completely different environmental conditions resulting in new challenges for the application of the corresponding methodologies. The need of high spatial sampling resolution for the different geophysical measurements to map efficiently the archaeological structures within a few meters below the seabed definitely leads to more effort and higher costs in relation to a respective land survey. However, this limitation is counterbalanced by the absence of an alternative geophysical method to extract valuable archaeological information in such an ultra-shallow aquatic environment.
The successful geophysical results from four different archaeological sites in the eastern Mediterranean can be used as a guide to compile a protocol or a sequence of procedures for the successful implementation of this methodological flowchart in diverse shallow marine environments worldwide. This approach will provide an extra management tool to the policy makers in order to plan and implement strategies for the documentation and conservation of shallow underwater Cultural Heritage. This is essentially crucial to preserve the history of humankind and ensure its accessibility to future generations.
Consequently, similar approaches could also raise the awareness for the local communities for more active participation. It is indeed strongly believed that only a proper and deep involvement of local community can assure a longer life to endangered submerged Cultural Heritage. Their attention and proudness for the proximity of the heritage will help in promoting informed tourism and support the development of the derived local economy in the longer term.
Funding: The compilation of this paper was supported by the project "DIATOPO-Cretan cultural landscapes over the time: highlighting the marine and mountainous environment of Mirabello" with MIS 5028190. DIATOPO is part of the Priority Axis "Strengthening the competitiveness, innovation and entrepreneurship of Crete" of the Operational Program, "Crete" 2014-2020 and co-financed by the European Regional Development Fund.
Institutional Review Board Statement: Ethical review and approval were waived for this study.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The spatial raw data collected in the archaeological sites are freely available to the academic and research community from the author upon reasonable request.