The Messapic Site of Muro Leccese: New Results from Integrated Geophysical and Archaeological Surveys

The regular application of geophysical survey techniques to evaluate archaeological sites is well established as a method for locating, defining, and mapping buried archaeological materials. However, it is not always feasible to apply a range of different methods over a particular site or landscape due to constraints in time or funding. This paper addresses the integrated application of three geophysical survey methods over an important archaeological site located in south Italy. In particular, it is focused on the results achieved from a past geophysical survey and the ongoing excavations performed by archaeologists in the site of Muro Leccese. Muro Leccese (Lecce) is one of the most important Messapian archaeological sites in southern Italy. The archaeological interest of the site was generated since the discovery of the remains of Messapian walls (late 4th–3rd centuries BC). With the aim of widening the archaeological knowledge of the Messapian settlement, several integrated methods, including magnetometry, ground-penetrating radar, and electrical resistivity tomography were used on site to fulfill a number of different research objectives. Since the most important targets were expected to be located at shallow soil depth, a three-dimensional (3D) ground-penetrating radar (GPR) survey was carried out in two zones, which were labeled respectively as zone 1 and zone 2, and were both quite close to the archaeological excavations. The GPR investigations were integrated with a 3D electrical resistivity tomography (ERT) survey in zone 1 and with a magnetometric, in gradiometry configuration survey in zone 2. The integration of several techniques allowed mapping the structural remains of this area and leading the excavation project. The geophysical results show a good correspondence with the archaeological features that were found after the excavation. Current work on the geophysical survey data using different codes for the processing of the data and merging different datasets using a Geographic Information System allowed achieving a user-friendly visualization that was presented to the archaeologists.

They are being increasingly utilized and reported upon by archaeological research teams for the detection and location of buried structures in terrestrial environments. In particular, several In particular, we have performed geophysical investigations in two zones labeled zone 1 and zone 2. The interpretation of GPR data were enhanced thanks to the ERT data in zone 1 and magnetic data in zone 2. The integrated results enabled us to provide a reasonable interpretation of the GPR data, which was also subsequently confirmed by some excavations, as it will be shown.

Site Description
The settlement of Muro Leccese, which is today a town of about 5000 inhabitants, is located in the Salento peninsula in the southern part of the Puglia Region, about 40 km southeast of Lecce ( Figure 1). Similar to other towns in this geographical area, its history began in the Iron Age (8th century BC), and has continued almost uninterrupted up to the present (Figure 1). The archaeological research conducted in Muro by the University of the Salento every year since 2000 has addressed several areas, albeit at different levels of detail, and the results now allow us to draw a preliminary profile of the ancient inhabited area as a whole [33].
Muro's long history can be divided into three phases characterized by differing kinds of settlement: During the Iron Age, from the 8th to the mid-6th centuries BC, the settlement was organized into clusters of huts distributed over an area of about 70 hectares. Inhabited by groups of interrelated families, the clusters were separated from each other by open spaces [34].
Remote Sens. 2019, 11, 1478 3 of 29 In particular, we have performed geophysical investigations in two zones labeled zone 1 and zone 2. The interpretation of GPR data were enhanced thanks to the ERT data in zone 1 and magnetic data in zone 2. The integrated results enabled us to provide a reasonable interpretation of the GPR data, which was also subsequently confirmed by some excavations, as it will be shown.

Site Description
The settlement of Muro Leccese, which is today a town of about 5000 inhabitants, is located in the Salento peninsula in the southern part of the Puglia Region, about 40 km southeast of Lecce ( Figure 1). Similar to other towns in this geographical area, its history began in the Iron Age (8 th century BC), and has continued almost uninterrupted up to the present (Figure 1). The archaeological research conducted in Muro by the University of the Salento every year since 2000 has addressed several areas, albeit at different levels of detail, and the results now allow us to draw a preliminary profile of the ancient inhabited area as a whole [33].
Muro's long history can be divided into three phases characterized by differing kinds of settlement: During the Iron Age, from the 8 th to the mid-6 th centuries BC, the settlement was organized into clusters of huts distributed over an area of about 70 hectares. Inhabited by groups of interrelated families, the clusters were separated from each other by open spaces [34]. In the course of the 6 th century BC (Archaic period), the town underwent profound transformations in terms of layout, architecture, and social structure, which were influenced by the nearby Greek colonies of southern Italy: the huts were abandoned in favor of houses made of stone and bricks, the banquet was adopted as a ceremonial practice reserved for the aristocracy, funerary rituals envisaged the burial of the corpse, and the Greek alphabet was used for inscriptions in the Messapian language [35][36][37].
At the end of the 4 th century BC, the landscape was marked by the construction of an imposing circuit of walls built with large squared blocks. Almost 4 km long and 3 meters thick, the walls In the course of the 6th century BC (Archaic period), the town underwent profound transformations in terms of layout, architecture, and social structure, which were influenced by the nearby Greek colonies of southern Italy: the huts were abandoned in favor of houses made of stone and bricks, the banquet was adopted as a ceremonial practice reserved for the aristocracy, funerary rituals envisaged the burial of the corpse, and the Greek alphabet was used for inscriptions in the Messapian language [35][36][37].
At the end of the 4th century BC, the landscape was marked by the construction of an imposing circuit of walls built with large squared blocks. Almost 4 km long and 3 m thick, the walls enclosed an area of more than 100 hectares [38]. The archaeological investigation revealed how its construction changed Remote Sens. 2019, 11,1478 4 of 26 the shape of the ancient settlement, the walls resetting the boundaries of the town. Some of the clusters of dwellings dating from the Archaic period lay outside the circuit and were abandoned [35,37,38].
However, the new walled settlement did not endure for a long time, and the town was almost certainly conquered and destroyed by the Romans following a siege. There are evidences of fires and collapses, followed by the rapid and definitive abandonment of the houses [39].
From a geological point of view, the Salento peninsula including the investigated archaeological area is the southernmost emerged part of Adria Plate, which constitutes the foreland of both Apenninic and Dinaric origins. It comprises a Variscan basement covered by a 3-5-km thick Mesozoic carbonate sequence (the "Calcari delle Murge" unit) overlain by thin Tertiary and Quaternary deposits. The most ancient rocks of this coverage were produced by transgressions after the definitive emersion of the Apulian carbonate platform, which occurred between the end of the Cretaceous and the beginning of Paleogene periods. In some cases, bauxitic deposits can be found between Mesozoic limestones and Paleogene units [19]. Four sedimentary cycles have been recognized from Neogene to Lower Pleistocene periods [3,19]. The first cycle comprises the formation of the "Pietra Leccese" and the overlying Calcarenite di Andrano Formation.
The Miocene sedimentary cycle was interrupted because of the emersion of Salento, which prevented the formation of Messinian evaporites. The total thickness of Miocene formations is greater than 150 m on the eastern side of the peninsula. The second cycle is represented by breccias and conglomerates of the Leuca Formation, which were deposited during the Lower Pliocene, reaching a maximum thickness of 30 m. The third sedimentary cycle is represented by the Upper Pliocene "Uggiano la Chiesa" Formation, which is composed of well-stratified and fossiliferous biodetritical limestones and yellowish calcareous sands with a maximum thickness of about 80 m. The fourth sedimentary cycle promoted the deposition in the Lower Pleistocene of the Calcareniti del Salento Formation, which is a very fossiliferous biodetritical calcareous sediment marked by the occurrence of Arctica islandica Linneo. Its maximum thickness is about 60 m. Finally, a number of Middle-Upper Pleistocene deposits related to eustatic sea level change can be found too on the Salento peninsula [3,19]. Three main tectonic events affected the Salento Peninsula during the Eo-Oligocene, the Middle Pliocene, and the Middle Pleistocene periods. In particular, the most recent tectonic phase was responsible for the final uplift of the Apulia foreland after the general subsidence took place in the early Miocene period; uplift ended at marine isotope stages (MIS) 9.3, about 330 ka BP. The Salento peninsula is a low-elevated landscape composed of a number of Pleistocene plains placed at different altitudes between sea level and 160 m of elevation. They are bordered by geomorphically degraded fault scarps, which were mostly elongated in the NW−SE and NNW−SSE directions, by differential erosion scarps, and by relict cliffs. Along the coast of the Salento peninsula, a number of marine terraces produced by the superimposition of regional uplift and glacio-eustatic sea level changes that occurred since the Middle Pleistocene period, can be recognized. In particular, the inner and western parts of the Salento Peninsula show a landscape produced by contact karst processes made of wide karstic surfaces, remnants of a Middle Pleistocene sedimentary cover, and morphostructural ridges. The karstic surfaces belong to a re-exhumed landscape shaped between the Lower and the Middle Pleistocene, whereas the morphostructural ridges are made of Mesozoic dolomitic-carbonatic units and show a polycyclic landscape. However, sinkholes are the main karstic landforms in many sectors of Salento, affecting all the outcropping carbonate rocks, including the Cretaceous limestone, the Oligocene, Miocene, and Plio-Pleistocene calcarenites, and the Middle and Upper Pleistocene terraced marine deposits. They are widespread, especially along the low elevated rocky platforms occurring on both the Adriatic and Ionian coasts, and display different typologies and states of activity. Locally, in the site at hand, data from archaeological excavations and geological observations (box core samples) made it possible to reconstruct the stratigraphy of the first 10 m of sediments. On the surface, there is a layer that is about 0.3 m thick constituted by covering materials (such as agricultural terrain). Beneath this layer, an organogenic calcarenites, compact detrital limestones layer is present, reaching up to 5 m in depth. layer is present, reaching up to 5 m in depth. In some places, this calcarenite is covered by up to three meters of alternating sands, silty sands, and marly limestone deposits (Figures 2 and 3) [32]. As can be seen from Figure 2, Muro Leccese is just on the boundary between two different geological formations, which made it difficult to predict a priori the possible characteristics of the soil. On the other hand, there was abundant archaeological information, as said previously. In particular, in 2002, the use of geophysical methods in the close Masseria Cunella district led to the detection of a possible grave burial. The following excavation campaign of 2003, conducted by Liliana Gardino et al., discovered a funerary area with two semi-chamber tombs that form a single complex used by the same family from the late 6 th to the mid-3 rd century B.C., with the deposition of 27 individuals. The grave goods highlight the occupants' high economic level, which in particular was witnessed by a late-Archaic Attic krater from the workshop of the Painter of Antimenes that was found in Tomb 2 [40][41][42]. As can be seen from Figure 2, Muro Leccese is just on the boundary between two different geological formations, which made it difficult to predict a priori the possible characteristics of the soil. layer is present, reaching up to 5 m in depth. In some places, this calcarenite is covered by up to three meters of alternating sands, silty sands, and marly limestone deposits (Figures 2 and 3) [32]. As can be seen from Figure 2, Muro Leccese is just on the boundary between two different geological formations, which made it difficult to predict a priori the possible characteristics of the soil. On the other hand, there was abundant archaeological information, as said previously. In particular, in 2002, the use of geophysical methods in the close Masseria Cunella district led to the detection of a possible grave burial. The following excavation campaign of 2003, conducted by Liliana Gardino et al., discovered a funerary area with two semi-chamber tombs that form a single complex used by the same family from the late 6 th to the mid-3 rd century B.C., with the deposition of 27 individuals. The grave goods highlight the occupants' high economic level, which in particular was witnessed by a late-Archaic Attic krater from the workshop of the Painter of Antimenes that was found in Tomb 2 [40][41][42]. On the other hand, there was abundant archaeological information, as said previously. In particular, in 2002, the use of geophysical methods in the close Masseria Cunella district led to the detection of a possible grave burial. The following excavation campaign of 2003, conducted by Liliana Gardino et al., discovered a funerary area with two semi-chamber tombs that form a single complex used by the same family from the late 6th to the mid-3rd century B.C., with the deposition of 27 individuals. The grave goods highlight the occupants' high economic level, which in particular was witnessed by a late-Archaic Attic krater from the workshop of the Painter of Antimenes that was found in Tomb 2 [40][41][42].
Following such an extraordinary result, other areas were investigated through geophysical methods before the archaeological excavation, namely Proprietà Natali (zone 1) and Palombara district (zone 2).

Zone 1
Zone 1 was divided in five areas labeled respectively A, B, C, D, and E ( Figure 4a). GPR propsecting was performed in all the areas, while ERT measurements were performed only in area C.
A GPR survey was performed using a pulsed RIS Hi Mod GPR system manufactured by the IDS-Corporation. The system is equipped with a dual-band antenna with central frequencies at 200 and 600 MHz, although only the results from the high-resolution 600 MHz survey are presented here. The archaeological targets looked for were expected to have a size between 0.5-2 m, and their depth was expected to range from tens of centimeters to a few meters. Following such an extraordinary result, other areas were investigated through geophysical methods before the archaeological excavation, namely Proprietà Natali (zone 1) and Palombara district (zone 2).

Zone 1
Zone 1 was divided in five areas labeled respectively A, B, C, D, and E ( Figure 4a). GPR propsecting was performed in all the areas, while ERT measurements were performed only in area C.
A GPR survey was performed using a pulsed RIS Hi Mod GPR system manufactured by the IDS-Corporation. The system is equipped with a dual-band antenna with central frequencies at 200 and 600 MHz, although only the results from the high-resolution 600 MHz survey are presented here. The archaeological targets looked for were expected to have a size between 0.5-2 m, and their depth was expected to range from tens of centimeters to a few meters. The survey was carried out in two zones (1 and 2 in Figure 4), which were both placed on the east side of the excavation. The GPR data were collected along parallel profiles spaced 0.5 m from each other. The GPR data were subsequently processed using standard two-dimensional processing techniques by means of the GPR-Slice Version 7.0 software [43]. The processing flow chart consisted of the following steps: (1) header editing for inserting the geometrical information; (2) frequency filtering; (3) manual gain to adjust the acquisition gain function and enhance the visibility of deeper anomalies; (4) customized background removal to attenuate the horizontal banding in the deeper part of the sections (ringing), performed by subtracting in different time ranges a local average trace estimated from suitably selected time-distance windows with low signal content (this local subtraction procedure was necessary to avoid artefacts created by the classic subtraction of the global average trace estimated from the entire section, due to the presence of zones with a very strong signal); (5) estimation of the average electromagnetic wave velocity by hyperbola fitting; and (6) Kirchhoff migration.
In Zone 1 area C, GPR measurements were integrated with ERT measurements. In particular, a special ERT arrays was used, with the electrodes distributed along a path surrounding area C [15,30,44] along an L-shaped line. Moreover, a dipole-dipole equatorial-parallel array was used, too. Initially, a 2D survey was conducted along each perpendicular line or transect, and subsequently, the current electrodes remained at the end of one line, while the potential electrodes were moved along the line. Then, the current electrodes were moved with a spatial step equal to the distance between The survey was carried out in two zones (1 and 2 in Figure 4), which were both placed on the east side of the excavation. The GPR data were collected along parallel profiles spaced 0.5 m from each other. The GPR data were subsequently processed using standard two-dimensional processing techniques by means of the GPR-Slice Version 7.0 software [43]. The processing flow chart consisted of the following steps: (1) header editing for inserting the geometrical information; (2) frequency filtering; (3) manual gain to adjust the acquisition gain function and enhance the visibility of deeper anomalies; (4) customized background removal to attenuate the horizontal banding in the deeper part of the sections (ringing), performed by subtracting in different time ranges a local average trace estimated from suitably selected time-distance windows with low signal content (this local subtraction procedure was necessary to avoid artefacts created by the classic subtraction of the global average trace estimated from the entire section, due to the presence of zones with a very strong signal); (5) estimation of the average electromagnetic wave velocity by hyperbola fitting; and (6) Kirchhoff migration.
In Zone 1 area C, GPR measurements were integrated with ERT measurements. In particular, a special ERT arrays was used, with the electrodes distributed along a path surrounding area C [15,30,44] along an L-shaped line. Moreover, a dipole-dipole equatorial-parallel array was used, too. Initially, a 2D survey was conducted along each perpendicular line or transect, and subsequently, the current electrodes remained at the end of one line, while the potential electrodes were moved along the line. Then, the current electrodes were moved with a spatial step equal to the distance between two adjacent electrodes, and the potential electrodes were again moved, as previously described. The process was repeated until the current and potential electrodes covered the L geometry. This sequence of observations produced a series of apparent resistivity observations toward and beneath the central portion of the array. The process is discussed in detail in [30].
The resistivity data were collected using a Syscal kid swich device (IRIS Instruments, France) supporting 24 electrodes. The electrode separation was 1 m. A penetration depth of 3.6 m was obtained. After the data acquisition, processing was performed, and the apparent resistivity data were analyzed to identify abnormal measurements with a high standard deviation.
The investigated volume was computed using the software ErtLab (http://www.geostudiastier.it), which makes use of finite element method (FEM). Several horizontal depth slices were analyzed to observe the extension versus depth of the resistivity features.

Zone 2
In the Zone 2, GPR and magnetic gradiometry data were acquired. GPR data were acquired and processed in the same way as described above for Zone 1. Magnetic gradiometry allows achieving a rapid mapping of magnetized archaeological objects, buried structures, and features [17]. Generally speaking, it highlights buried structures that show an appreciable magnetic susceptibility and therefore high contrast with the magnetic susceptibility of the embedding soils and sediments. Standard gradiometers measure the vertical gradient of the local magnetic field. The values of the vertical gradient are influenced by shallow magnetic subsurface archaeological features. The depth of investigation using the gradiometry technique depends on the nature of possible subsurface targets and can extend up to 2 m.
The objective of the gradiometric survey is to detect, by means of surface measurements, variations in the magnetic properties, or magnetization, of the subsurface. These variations in subsurface magnetization can be related to archaeological remains (walls, tombs, floors, etc.) as well as geological structures such as a pegmatite vein or basaltic dike in granite, a fault in bedrock, or hydrothermal alterations. They can also be related to the human disturbance of soil through the presence of magnetic objects in the subsurface. The anomalous anomalies define the spatial variations in the field that can be measured on the surface. Mapping magnetic anomalies on the surface allows inferring the presence or absence of magnetic material in the subsurface. The magnetic gradient allows differentiation between deeply buried objects and shallower ones [17]. Typically, a gradiometer uses two sensors to measure the gradient in the magnetic field. In the case history at hand, the data were acquired using a Bartington model 601 gradiometer along parallel lines spaced 1 m apart. Data were acquired in continuous mode with a sampling interval of 0.2 s. The sensors were oriented vertically. A data logger and a control console, which was mounted in front of the driver, were used to control the data collection. Data were recorded with an accuracy of 0.002 nano-Teslas (nT). All the acquired magnetic data were preprocessed with the use of bandpass amplitude filters with values between −30nT/m+30nT/m, which represents the range wherein the values associated to archaeologically significant magnetic alignments usually fall. In this phase, in addition, it was necessary to eliminate events, which are referred to as spikes, and were attributable to the presence of alignments due to anthropogenic settlement. In the next step, all the data were interpolated through a grid of data scaled between −20 nT/m and +20 nT/m.  Some strong reflection events (labeled M) have been interpreted as due to walls, while other reflection events (labeled P) have been interpreted as due to pipes ( Figure 5).

Results
The processed GPR data were subsequently merged together into three-dimensional volumes and visualized in various ways in order to enhance the spatial correlations among the anomalies of interest. A way to obtain useful maps for understanding the plan distribution of reflection amplitudes within specific time intervals is the creation of horizontal time slices. These are maps on which the reflection amplitudes are projected at a specified time (or depth), averaged within a selected time interval [31]. In a graphic method developed by [45] termed ''overlay analysis'', the strongest and weakest reflectors within each slice are enhanced with specific colors. This technique allows the linkage of structures buried at different depths. This represents an improvement in imaging, because subtle features that are indistinguishable in the radargrams in this way can be more easily seen and interpreted. In the present work, the time-slice technique has been used to display the amplitude variations within consecutive time windows of width Δt = 5 ns. Strong reflection events (labeled P and M) are clearly visible ( Figure 6). Some strong reflection events (labeled M) have been interpreted as due to walls, while other reflection events (labeled P) have been interpreted as due to pipes ( Figure 5).
The processed GPR data were subsequently merged together into three-dimensional volumes and visualized in various ways in order to enhance the spatial correlations among the anomalies of interest. A way to obtain useful maps for understanding the plan distribution of reflection amplitudes within specific time intervals is the creation of horizontal time slices. These are maps on which the reflection amplitudes are projected at a specified time (or depth), averaged within a selected time interval [31]. In a graphic method developed by [45] termed "overlay analysis", the strongest and weakest reflectors within each slice are enhanced with specific colors. This technique allows the linkage of structures buried at different depths. This represents an improvement in imaging, because subtle features that are indistinguishable in the radargrams in this way can be more easily seen and interpreted. In the present work, the time-slice technique has been used to display the amplitude variations within consecutive time windows of width ∆t = 5 ns. Strong reflection events (labeled P and M) are clearly visible ( Figure 6). A further approach for visualizing 3D images is the use of surface rendering, assembled from all the processed radargrams [18,20,27].
In Figure 7, the same GPR data set is displayed with isoamplitude surfaces using a threshold value of 70% for the maximum complex trace amplitude. Whilst lowering the threshold value increases the visibility of the main anomaly and smaller objects, it also increases the undesired effects A further approach for visualizing 3D images is the use of surface rendering, assembled from all the processed radargrams [18,20,27].
In Figure 7, the same GPR data set is displayed with isoamplitude surfaces using a threshold value of 70% for the maximum complex trace amplitude. Whilst lowering the threshold value increases the visibility of the main anomaly and smaller objects, it also increases the undesired effects due to the present heterogeneities. Relatively strong continuous reflection events are visible on the thresholded volumes (P and M). This visualization technique provides strong evidence for the anomaly related to the pipe (P) and wall (M). A further approach for visualizing 3D images is the use of surface rendering, assembled from all the processed radargrams [18,20,27].
In Figure 7, the same GPR data set is displayed with isoamplitude surfaces using a threshold value of 70% for the maximum complex trace amplitude. Whilst lowering the threshold value increases the visibility of the main anomaly and smaller objects, it also increases the undesired effects due to the present heterogeneities. Relatively strong continuous reflection events are visible on the thresholded volumes (P and M). This visualization technique provides strong evidence for the anomaly related to the pipe (P) and wall (M).  The reflection event shown in Figure 8b appears to be of particular interest. In fact, there is an evident polarity inversion of the reflected electromagnetic wave. This event (T) could be related to a strong contrast of electromagnetic characteristic due for example to a presence of a void room [20,21,28]. Therefore, this event was interpreted as due to the probable presence of a tomb. The reflection event shown in Figure 8b appears to be of particular interest. In fact, there is an evident polarity inversion of the reflected electromagnetic wave. This event (T) could be related to a strong contrast of electromagnetic characteristic due for example to a presence of a void room [20,21,28]. Therefore, this event was interpreted as due to the probable presence of a tomb.
The time-slice technique allowed visualizing the extension of the reflected events labeled A, M, and T, respectively, in Figure 9. The isoamplitude 3D reconstruction achieved using a threshold value of 70% for the maximum complex trace amplitude is shown in Figure 10. A relatively strong continuous reflection event is visible on the thresholded volumes (see anomalies A, M, and T).
.  The isoamplitude 3D reconstruction achieved using a threshold value of 70% for the maximum complex trace amplitude is shown in Figure 10. A relatively strong continuous reflection event is visible on the thresholded volumes (see anomalies A, M, and T). The isoamplitude 3D reconstruction achieved using a threshold value of 70% for the maximum complex trace amplitude is shown in Figure 10. A relatively strong continuous reflection event is visible on the thresholded volumes (see anomalies A, M, and T).
.  The hyperbola analysis showed also an increasing of the electromagnetic wave propagation velocity from 0.08 m/ns (esteemed in the areas A and B) to 0.112 m/ns. In this area, a particular interest seems to have the reflection event shown in Figure 11a,b. A polarity inversion of the reflected electromagnetic wave is also evident in this area, which is in relationship with the event labeled A (Figure 11). The time-slice technique allows visualizing the extension of the reflected event labeled as A in Figure 12. The isoamplitude 3D perspective reconstruction using a threshold value 70% for the maximum complex trace amplitude that is shown in Figure 13 puts into evidence the dimension of the event labeled A. Figure 11. Zone 1, area C: an example of the two processed radar sections (a,b). It is possible to see "A" large tomb This event (A) was interpreted as due to a probable presence of large tomb (about 4 × 2.5 m 2 ). The time-slice technique allows visualizing the extension of the reflected event labeled as A in Figure 12. The isoamplitude 3D perspective reconstruction using a threshold value 70% for the maximum complex trace amplitude that is shown in Figure 13 puts into evidence the dimension of the event labeled A.    In order to properly understand the nature of the anomaly labeled "A" in the GPR results, ERT measurements were also performed in area C. Figure 14 shows ERT slices from 0.5 to 3.0 m in depth. In order to properly understand the nature of the anomaly labeled "A" in the GPR results, ERT measurements were also performed in area C. Figure 14 shows ERT slices from 0.5 to 3.0 m in depth. It is possible to note the presence of a heterogeneous subsurface with resistivity values ranging from 2000 to 4000 Ω m −1 . Then, it is also possible to note the presence of the area indicated with ''A'', with resistivity values between 3500-4000 Ω m −1 ; these values indicate the probable presence of a void. Figure 15 shows the 3D contouring of isoresistivity surfaces. In this representation, the transparency function is defined by two threshold values for the resistivity ρ1 and ρ2 (ρ1 < ρ2). In the intervals ρ < ρ1 and ρ > ρ2, the data is made transparent; therefore, only the data in the interval ρ1 < ρ < ρ2 are visualized. The threshold calibration is a very delicate task. In fact, by lowering the threshold value, not only is the visibility of the main anomaly raised, but that of the smaller objects and noise also increases. In Figure 15, the resistivity data set is displayed with isoresistivity surfaces using a threshold value ranging from 3000 to 4000 Ω m. In Figure 15b, the resistivity data set is displayed using the horizontal and vertical slices. The lowest resistivity anomalies were removed, and therefore, the cavity (A) spaces are well visible. It is possible to note the presence of a heterogeneous subsurface with resistivity values ranging from 2000 to 4000 Ω m −1 . Then, it is also possible to note the presence of the area indicated with "A", with resistivity values between 3500-4000 Ω m −1 ; these values indicate the probable presence of a void. Figure 15 shows the 3D contouring of isoresistivity surfaces. In this representation, the transparency function is defined by two threshold values for the resistivity ρ1 and ρ2 (ρ1 < ρ2). In the intervals ρ < ρ1 and ρ > ρ2, the data is made transparent; therefore, only the data in the interval ρ1 < ρ < ρ2 are visualized. The threshold calibration is a very delicate task. In fact, by lowering the threshold value, not only is the visibility of the main anomaly raised, but that of the smaller objects and noise also increases. In Figure 15, the resistivity data set is displayed with isoresistivity surfaces using a threshold value ranging from 3000 to 4000 Ω m. In Figure 15b, the resistivity data set is displayed using the horizontal and vertical slices. The lowest resistivity anomalies were removed, and therefore, the cavity (A) spaces are well visible.  It is evident that a classical metal reflection labeled "P" related to a pipe. The reflection events (M) were interpreted as due to the probable presence of walls. In particular, the reflection event M1 seem to be due to pebbles put together to form a wall. The time-slice technique allows visualizing the extension of the reflected events labeled P, M, and M1 ( Figure 17). It is evident that a classical metal reflection labeled "P" related to a pipe. The reflection events (M) were interpreted as due to the probable presence of walls. In particular, the reflection event M1 seem to be due to pebbles put together to form a wall. The time-slice technique allows visualizing the extension of the reflected events labeled P, M, and M1 ( Figure 17). It is evident that a classical metal reflection labeled "P" related to a pipe. The reflection events (M) were interpreted as due to the probable presence of walls. In particular, the reflection event M1 seem to be due to pebbles put together to form a wall. The time-slice technique allows visualizing the extension of the reflected events labeled P, M, and M1 ( Figure 17).  The isoamplitude surfaces using the threshold value of 67% for the maximum complex trace amplitude shown in Figure 18 puts into evidence the dimensions of the events labeled P, M, and M1.  (Figure 19a,b). It is evident that a "cave-like" reflection labeled "T" related to a possible tomb. The reflection events (M) were interpreted as due to a probable presence of walls. The time-slice technique allows visualizing the extension of the reflected events labeled T and M ( Figure 20). The reflection event "W" was interpreted as wall foundations.  (Figure 19a,b). It is evident that a "cave-like" reflection labeled "T" related to a possible tomb. The reflection events (M) were interpreted as due to a probable presence of walls. The time-slice technique allows visualizing the extension of the reflected events labeled T and M ( Figure 20). The reflection event "W" was interpreted as wall foundations.

Zone 1, Area E
Area E was sized about 20 × 22 m. The hyperbola analysis showed an electromagnetic wave propagation velocity of 0.094 m/ns. In this area, several reflection events are visible (Figure 19a,b). It is evident that a "cave-like" reflection labeled "T" related to a possible tomb. The reflection events (M) were interpreted as due to a probable presence of walls. The time-slice technique allows visualizing the extension of the reflected events labeled T and M (Figure 20). The reflection event "W" was interpreted as wall foundations.  The isoamplitude surface 3D reconstruction using the threshold value of 65% for the maximum complex trace amplitude shown in Figure 21 puts into evidence the dimensions of the events labeled W, M, and T. The isoamplitude surface 3D reconstruction using the threshold value of 65% for the maximum complex trace amplitude shown in Figure 21 puts into evidence the dimensions of the events labeled W, M, and T. The isoamplitude surface 3D reconstruction using the threshold value of 65% for the maximum complex trace amplitude shown in Figure 21 puts into evidence the dimensions of the events labeled W, M, and T.

Zone 2
Zone 2 was a square sized 40 × 40 m and was divided in four square areas, each of which were sized 20 × 20 m, and were respectively labeled F, G, H, and I ( Figure 4). Ground-penetrating radar (GPR) and magnetometric geophysical measurements in gradiometry configuration were used. In Zone 2, the gradiometry magnetic method was used almost everywhere, while the GPR method was used only in area F.

Zone 2
Zone 2 was a square sized 40 × 40 m and was divided in four square areas, each of which were sized 20 × 20 m, and were respectively labeled F, G, H, and I ( Figure 4). Ground-penetrating radar (GPR) and magnetometric geophysical measurements in gradiometry configuration were used. In Zone 2, the gradiometry magnetic method was used almost everywhere, while the GPR method was used only in area F.  As said, area F was sized 20 × 20 m 2 . A Ris Hi-mod georadar system with a dual band 600 MHz-200 MHz antennas was used. The GPR data were acquired with a transect of 0.5 m between any two adjacent measurement lines, and the data were processed by making use of the GPRSlice code. In

Zone 2: Area F, GPR Data Analysis
As said, area F was sized 20 × 20 m 2 . A Ris Hi-mod georadar system with a dual band 600 MHz-200 MHz antennas was used. The GPR data were acquired with a transect of 0.5 m between any two adjacent measurement lines, and the data were processed by making use of the GPRSlice code. In particular, here, we focus on the results obtained with the 600-MHz antenna. With reference to Figure 23, the analysis of the GPR data pointed out the presence of several anomalies ("M") that represent walls; Anomaly C, instead, indicates the probable presence of a tomb.  Slices analysis (Figure 24) show the distribution of anomalies "M" and "C". The 3D isosurfaces amplitude representation, with a threshold of 70% (Figure 25a,b), show the 3D distribution of the "M" and "C" anomalies.
It might have to be noted at this point that the threshold for the isosurface representation was not the same throughout the paper. Indeed, an optimal value cannot be established a priori in general, and the choice is necessarily case-dependent. It is a good role of thumb to check that the main anomalies visible in the isosurface representation are robust with respect to the choice of the threshold [46]. We have done it, even if for sake of brevity we have avoided presenting all of the synopses here. The 3D isosurfaces amplitude representation, with a threshold of 70% (Figure 25a,b), show the 3D distribution of the "M" and "C" anomalies.  The 3D isosurfaces amplitude representation, with a threshold of 70% (Figure 25a,b), show the 3D distribution of the "M" and "C" anomalies.

Archaeological Excavations in Zone 1
In July 2012, two archaeological excavations were made in the Natali estate, within zone 1 (Figures 1 and 2): the first one was in Area B, and the second one was in Area E. Since this is a private property, these excavations were conducted in the areas of the field not cultivated or built up: in area C, instead, the presence of a base of cement made by the owner did not allow digging and verifying the results of the geophysical survey.

Excavation in Area B
The 5 × 3 m archaeological excavation in this area ( Figure 26) has revealed the presence, in its eastern part, of the foundation of a wall made up of roughly squared sandstone blocks of considerable size, oriented NE-SW. The blocks have evident mole plough traces on their upper face. The type of stone and the workmanship [47], as well as the retrieved archaeological materials, allow us dating it to the Hellenistic period, in the second half of the 4th century BC.
It might have to be noted at this point that the threshold for the isosurface representation was not the same throughout the paper. Indeed, an optimal value cannot be established a priori in general, and the choice is necessarily case-dependent. It is a good role of thumb to check that the main anomalies visible in the isosurface representation are robust with respect to the choice of the threshold [46]. We have done it, even if for sake of brevity we have avoided presenting all of the synopses here.

Archaeological Excavations in Zone 1
In July 2012, two archaeological excavations were made in the Natali estate, within zone 1 (figures 1 and 2): the first one was in Area B, and the second one was in Area E. Since this is a private property, these excavations were conducted in the areas of the field not cultivated or built up: in area C, instead, the presence of a base of cement made by the owner did not allow digging and verifying the results of the geophysical survey.

Excavation in Area B
The 5 × 3 m archaeological excavation in this area ( Figure 26) has revealed the presence, in its eastern part, of the foundation of a wall made up of roughly squared sandstone blocks of considerable size, oriented NE-SW. The blocks have evident mole plough traces on their upper face. The type of stone and the workmanship [47], as well as the retrieved archaeological materials, allow us dating it to the Hellenistic period, in the second half of the 4 th century BC. Immediately to the west of this structure, there is a thick pavement of compacted limestone dust that overlaps a previous cobblestone. The materials found in the layer between the two floors enable us to date the cobblestone to the Iron Age and relate it to the presence of a hut. To the east of the wall, instead, a layer of a reddish brown soil is present with materials dating back to the Iron Age (preliminary data about the Iron age pottery and structures are in [38].
Almost behind the wall, the Hellenistic floor has been cut by a semi-circular hole due to a tree, which was removed by the owner a few years before the excavation. Immediately to the west of this structure, there is a thick pavement of compacted limestone dust that overlaps a previous cobblestone. The materials found in the layer between the two floors enable us to date the cobblestone to the Iron Age and relate it to the presence of a hut. To the east of the wall, instead, a layer of a reddish brown soil is present with materials dating back to the Iron Age (preliminary data about the Iron age pottery and structures are in [38].
Almost behind the wall, the Hellenistic floor has been cut by a semi-circular hole due to a tree, which was removed by the owner a few years before the excavation.
Collecting these data, it is possible to affirm that the excavation has enabled highlighting a room of a Hellenistic house (second half of the 4th-3rd century BC), below which there was another dwelling, dating back to the Iron Age (8th-half of the 6th century BC).

Excavation in Area E
A surface of 7.50 × 4 m was excavated in this area (Figure 27a,b). Below the top layer of agricultural soil, which was on average 40-45 cm deep, the western edge of a road, made of small stones compacted with limestone dust, emerged in the eastern part of the excavation.
Remote Sens. 2019, 11, 1478 23 of 29 Collecting these data, it is possible to affirm that the excavation has enabled highlighting a room of a Hellenistic house (second half of the 4 th -3 rd century BC), below which there was another dwelling, dating back to the Iron Age (8 th -half of the 6 th century BC).

Excavation in Area E
A surface of 7.50 × 4 m was excavated in this area (Figure 27a,b). Below the top layer of agricultural soil, which was on average 40-45 cm deep, the western edge of a road, made of small stones compacted with limestone dust, emerged in the eastern part of the excavation. To the west of the road edge, a ramp-sidewalk about 2.5 m wide in well-purified white compacted limestone dust that links the road with a long wall structure has been brought to light. It was probably part of a dwelling, and in particular, the eastern portion of two rooms.
Between the edge of the road and the edge of the ramp, there is a water channel that is approximately 25 cm wide for draining rainwater. Even in the ramp itself, there is a long transverse cut, which was probably linked to another channel for the drainage of water from the building toward the road.
The two rooms of the probable dwelling have been excavated only in their northern edge. The type of stone of the blocks and the construction technique allow us to date the structure to the Hellenistic period. Under a thin slab of limestone dust floor, in a poor state of preservation, a level with material dating back from the Iron Age to the archaic period was found. A layer of natural soil or the non-regular bedrock lies under it, depending on the points.
The wall structures discovered during the excavation of areas B and E and the possible wall identified by the geophysical investigations in area A (figures 5 and 7) are perpendicular and parallel to each other. It suggests that they can belong to a large residential dwelling, which was most probably similar to those found in the Cunella [33,39,41,[47][48][49] and Palombara [38] districts, with the front of the building aligned along the road that led toward the north gate of the surrounding walls.

Archaeological Excavation in Zone 2
During the excavation periods of July and September 2016 in the Palombara district (zone 2), a surface of 3 × 5 m was excavated in area F to verify an anomaly discovered during the geophysical analyses and interpreted as a possible burial. The excavation was carried out in parallel with the To the west of the road edge, a ramp-sidewalk about 2.5 m wide in well-purified white compacted limestone dust that links the road with a long wall structure has been brought to light. It was probably part of a dwelling, and in particular, the eastern portion of two rooms.
Between the edge of the road and the edge of the ramp, there is a water channel that is approximately 25 cm wide for draining rainwater. Even in the ramp itself, there is a long transverse cut, which was probably linked to another channel for the drainage of water from the building toward the road.
The two rooms of the probable dwelling have been excavated only in their northern edge. The type of stone of the blocks and the construction technique allow us to date the structure to the Hellenistic period. Under a thin slab of limestone dust floor, in a poor state of preservation, a level with material dating back from the Iron Age to the archaic period was found. A layer of natural soil or the non-regular bedrock lies under it, depending on the points.
The wall structures discovered during the excavation of areas B and E and the possible wall identified by the geophysical investigations in area A (Figures 5 and 7) are perpendicular and parallel to each other. It suggests that they can belong to a large residential dwelling, which was most probably similar to those found in the Cunella [33,39,41,[47][48][49] and Palombara [38] districts, with the front of the building aligned along the road that led toward the north gate of the surrounding walls.

Archaeological Excavation in Zone 2
During the excavation periods of July and September 2016 in the Palombara district (zone 2), a surface of 3 × 5 m was excavated in area F to verify an anomaly discovered during the geophysical analyses and interpreted as a possible burial. The excavation was carried out in parallel with the excavation of the so-called inner city walls, which allowed highlighting the important remains of a late-Archaic and Classical age residential building aligned to a road five meters wide [36,38,50].
The removal of the top soil layer (Figure 28a,b), which was only a few centimeters deep, highlighted the bedrock, which was not revealed only in the area where the anomaly had been identified with the geophysical survey. excavation of the so-called inner city walls, which allowed highlighting the important remains of a late-Archaic and Classical age residential building aligned to a road five meters wide [36,38,50]. The removal of the top soil layer (Figure 28a,b), which was only a few centimeters deep, highlighted the bedrock, which was not revealed only in the area where the anomaly had been identified with the geophysical survey. A deep investigation at that point made it possible to verify that, beneath the top soil, there was an area of a whitish natural clay. Its excavation in depth showed the total absence of anthropic elements, and therefore the excavation was stopped at a depth of 2.6 m below the bedrock level.

Discussion and Conclusions
The research approach described and explained in this paper demonstrate how traditional archaeological studies can be enhanced through the application of multidisciplinary investigations.
The combination of GPR, ERT, and magnetic data allow a comprehensive interpretation of the subsurface structures present in the investigated site. The GPR was efficient within the upper 2 m to identify the different near-surface structures. In cases where the bedrock is at shallow depth, it is possible to distinguish it from unconsolidated materials on the GPR profiles.
The GPR investigations are based on the emission of electromagnetic waves, which propagate inside the ground and are scattered in all directions, including that coming back from the GPR when they meet changes in the dielectric properties of the materials constituting the subsoil. Generally, in archaeological applications, GPR methodology performs quite satisfyingly when highly resistive materials occur [20]. In fine-grained soils such as the investigated site, the presence of conductive materials is one of the most important factors determining GPR results. In fact, these materials impose a strong reduction of EM wave penetration inside the medium.
The ERT was suitable to achieve information about deeper buried structures. ERT is a geophysical technique for searching and imaging subsurface structures based on electrical resistivity measurements made at the surface using sets of metal electrodes.
The electrical properties of the subsurface vary with the ground material, the presence and saturation level of fluids, and the presence of buried objects. Electrical techniques retrieve the distribution of these properties as a function of depth and the horizontal distance [20]. Also for ERT, the presence of conductive materials is one of the most important factors determining the results. In fact, the conductive materials favor the passage of electric current and allow obtaining more detailed information about the buried structures. A deep investigation at that point made it possible to verify that, beneath the top soil, there was an area of a whitish natural clay. Its excavation in depth showed the total absence of anthropic elements, and therefore the excavation was stopped at a depth of 2.6 m below the bedrock level.

Discussion and Conclusions
The research approach described and explained in this paper demonstrate how traditional archaeological studies can be enhanced through the application of multidisciplinary investigations.
The combination of GPR, ERT, and magnetic data allow a comprehensive interpretation of the subsurface structures present in the investigated site. The GPR was efficient within the upper 2 m to identify the different near-surface structures. In cases where the bedrock is at shallow depth, it is possible to distinguish it from unconsolidated materials on the GPR profiles.
The GPR investigations are based on the emission of electromagnetic waves, which propagate inside the ground and are scattered in all directions, including that coming back from the GPR when they meet changes in the dielectric properties of the materials constituting the subsoil. Generally, in archaeological applications, GPR methodology performs quite satisfyingly when highly resistive materials occur [20]. In fine-grained soils such as the investigated site, the presence of conductive materials is one of the most important factors determining GPR results. In fact, these materials impose a strong reduction of EM wave penetration inside the medium.
The ERT was suitable to achieve information about deeper buried structures. ERT is a geophysical technique for searching and imaging subsurface structures based on electrical resistivity measurements made at the surface using sets of metal electrodes.
The electrical properties of the subsurface vary with the ground material, the presence and saturation level of fluids, and the presence of buried objects. Electrical techniques retrieve the distribution of these properties as a function of depth and the horizontal distance [20]. Also for ERT, the presence of conductive materials is one of the most important factors determining the results. In fact, the conductive materials favor the passage of electric current and allow obtaining more detailed information about the buried structures.
The magnetic method in a gradiometer configuration offers a high degree of immunity from diurnal and minor magnetic storm activity in the ambient magnetic field; they can enhance near-surface, small, or weakly magnetic anomalies. They can also provide improvements in spatial resolution over the total field measurement alone [19]. It allows locating small, near surface anomalies, and is therefore very useful in archeological mapping. It can be used in order to investigate wide areas when the expected results are the buried tombs (especially if they contain metal materials as e.g., swords, shell helmets, or also terra cotta materials or kilns).
For all of the GPR data collected in this study, those made with the 600-MHz antenna offered the best resolution, although attenuation caused lower penetration. On these profiles, it was possible to clearly identify the buried infrastructures in the area. In particular, in zone 1, a series of structures were evidenced in the first meter of depth. Deeper attention was given to the anomaly indicated with "A" in zone 1, area C. The electromagnetic wave reflected event, due to its shape and size, was interpreted as a cut in the rocky basement linked to the probable presence of a chamber tomb. However, it could be a simple cut in the rock. These difficulties in interpreting the GPR data have led to performing a geoelectrical survey in the same area. Therefore, the ERT results were introduced as a complement of GPR results. An area of high resistivity (>4500 Ω.m) was found to correspond to the anomaly labeled "A" in the GPR results. The resistivity survey was useful to determine the nature of the anomaly and confirmed the probable presence of the chamber tomb.
In zone 2, a magnetic survey was conducted in a gradiometric configuration. The magnetic gradiometry survey results were relatively good, evidencing some buried alignments. Some problems in the data interpretation were related to the depth of the buried structures. Only GPR data could definitively determine the depth of buried structures below ground level.
Finally, an upgrading and increment of the knowledge of the Messapian settlement of Muro Leccese has been obtained by the acquisition of new archaeological and geophysical data, which have been integrated in a new digital archaeological map (Figures 29 and 30). The research and the following excavations have documented the extension of the inhabited area, have attested the presence of a new residential building, and have increased the knowledge of the road map of the town. With regard to the necropoleis, the contribution of geophysical data has been very important, as it showed the possible presence of grave tombs, which will be investigated in the next few years.
All this work has been conducted in the framework of a proficuous interdisciplinary interaction between geophysical and archaeological competences, conducted before the prospecting, because the archaeologists have the sensibility about the questions that need to be addressed in the archaeological site, after which the geophysicists conducted the prospecting together with the archaeologists. Afterwards, the interpretations were discussed together, and the archaeological excavation performed and that still scheduled for the future have been planned to account for the results of the geophysical prospecting.