Automated Resistivity Profiling (ARP) to Explore Wide Archaeological Areas: The Prehistoric Site of Mont'e Prama, Sardinia, Italy

This paper deals with the resistivity continuous surveys on extensive area carried out at the Mont’e Prama archaeological site, in Sardinia (Italy). From 2013 to 2015, new research was performed using both non-destructive surveys and traditional archaeological excavations. The measurements were done in order to find geophysical anomalies related to unseen buried archaeological remains and to define the spatial extension of the ancient necropolis. The electrical resistivity of soils was measured by means of the Automated Resistivity Profiling (ARP©) system. This multi-pole method provided high-resolution maps of electrical resistivity in the whole investigated area using a computer-assisted acquisition tool, towed by a small vehicle. Through this acquisition layout, a surface of 22,800 m2 was covered. The electrical resistivity data were derived in real time with centimetric horizontal precision through a differential GPS positioning system. Thanks to the simultaneous acquisition of ARP and GPS data, the rigorous georeferencing of the tridimensional experimental dataset was made possible, as well as the reconstruction of a detailed Digital Terrain Model. Here, the experimental results are analyzed and critically discussed by means of the integration of the results obtained by a high-resolution prospection performed with a multi-channel Ground Penetrating Radar system and taking into account other information derived from previous geological and archaeological studies. Geophysical results, jointly with topographic reconstruction, clearly permitted the identification of more interesting areas where future archaeological investigations could be focused.


Introduction
Modern archaeological surveys require non-invasive techniques that can be employed on wide areas in order to study the shallow subsoil, highlighting the presence of physical anomalies that could be related to buried archaeological remains. These techniques should be able to identify the physical properties of the soils by means of measuring systems positioned on a mobile vehicle in order to reduce the time of acquisition compared to other more traditional methods [1].
Electrical methods based on fixed arrays permit the creation of resistivity models of the medium through a set of electrodes arranged along linear profiles or along two-dimensional configurations [23]. 2 of 22 In the case of two-dimensional Electrical Resistivity Tomographies (ERT) the surveys provide a vertical section of the medium. The distance between the electrodes is constant. Each measurement is related to a specific distance between current and potential electrodes. Loops of electrodes with various geometries are used to perform tri-dimensional electrical resistivity tomographies in order to obtain a 3D model of the investigated medium [24][25][26][27][28][29][30][31][32][33][34].
Different systems for measuring and mapping electrical resistivity have been developed in order to overcome the main limitations of the fixed-array systems [35]. In particular, multi-polar mobile systems allow the performance of continuous resistivity profiling without a fixed device with an effective reduction in the time of acquisition. The first solution proposed a quadrupolar mobile device, provided with wheels equipped with electrodes for measuring the electrical resistivity of the soil at a constant theoretical depth [36]. The depth of investigation increased by varying the space between the wheels using an extendable support. Different solutions based on multi-polar mobile systems were developed and described in Panissod et al. [37]. These authors proposed an instrument based on eight electrodes named MUlti-pole Continuous Electrical Profiling (MUCEP), positioned along a V-shaped array (also called "Vol-de-canards" array) and located on a mobile measurement system. This equipment was used to explore shallow layers up to 3 meters in depth. In the same context, the comparison among resistivity prospections, obtained by the use of electrodes arrays towed by a vehicle, showed both advantages and limitations based on Pole-Pole arrays with a single depth of investigation with regards to the multi-pole systems. Multi-polar configuration provides a three-dimensional model of the subsoil with a significant reduction in acquisition time [38]. Further studies [35] have improved these multi-pole systems in order to reduce the noise effects on the raw data. In particular, the authors propose the equipment known as Automated Resistivity Profiling (ARP). This acquisition instrument is composed of a transversal current dipole and of three transversal receiving dipoles, set at an increasing distance from the transmitting one. This electrode array is designed to obtain simultaneously measurements of apparent resistivity simultaneously at three different subsoil depths. Papadopoulos et al. [39,40] described a three-dimensional inversion of resistivity data acquired with the ARP system in order to obtain a reliable numerical model capable of describing the spatial distribution of the electrical resistivity in the subsoil of the investigated area. The inversion process described by the authors is based on the least squares method and establishes smoothness constraints in order to consider the instability of the model and the non-uniqueness of solutions [41].
Several studies described multi-pole acquisition systems with the application of precision agriculture and soil characterization [42][43][44][45], archaeological prospection [46], and non-destructive evaluation of pavement and concrete infrastructures [47]. The ARP system was used by some authors [48] to evaluate the reduction of soil characterization costs for viticulture: the authors showed how the ARP system represents an effective method to reduce the cost of soil analysis, collecting a minor number of soil samples for direct tests. Many results of archaeological surveys carried out with the ARP system are extensively described in [1].
Integrated prospection with complementary geophysical methods permits to better define the physical model of the investigated subsoil [49]. One of the most integrated approaches in archaeology is the combination of electrical and Ground Penetrating Radar (GPR) methods [19,50].
In this paper, we present the results of surveys carried out by means of the ARP method at the archaeological site of Mont'e Prama (Sardinia, Italy). Geophysical investigations were designed to support future excavation activities and to define the spatial extension of the site in the northern side of an area in which archaeological digs took place in the 1970s. The analysis of electrical resistivity anomalies, both due to the buried archaeological remains and the geomorphological and geo-pedological spatial variability, permitted the evaluation of the effectiveness of this method for the site of study. Resistivity data were compared with the results achieved by extensive high-resolution surveys performed with multi-channel GPR at the same site. The archaeological interpretation of the experimental data was supported by the integration of different complementary geophysical methods and GPS data.

The Archaeological Site of Mont'e Prama
The study area of Mont'e Prama is one of the most important archaeological areas of the western Mediterranean [51]. The site is located in the Sinis peninsula (western-central Sardinia, Italy) as shown in Figure 1. The first archaeological remains were accidentally discovered in 1974 during soils ploughing when a local farmer found a fragment of a Nuragic sculpture. In 1975In , 1977In and 1979 various campaigns of archaeological excavation were performed and a necropolis with characteristic monumental sculptures from the early Iron Age was studied. Many baetyls, models of "nuraghe" and numerous fragments of sculptures (usually called "Giants") were collected [51,52]. The last archaeological trials of these campaigns did not find new evidence and excavations were suspended.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 23 experimental data was supported by the integration of different complementary geophysical methods and GPS data.

The Archaeological Site of Mont'e Prama
The study area of Mont'e Prama is one of the most important archaeological areas of the western Mediterranean [51]. The site is located in the Sinis peninsula (western-central Sardinia, Italy) as shown in Figure 1. The first archaeological remains were accidentally discovered in 1974 during soils ploughing when a local farmer found a fragment of a Nuragic sculpture. In 1975In , 1977 and 1979 various campaigns of archaeological excavation were performed and a necropolis with characteristic monumental sculptures from the early Iron Age was studied. Many baetyls, models of "nuraghe" and numerous fragments of sculptures (usually called "Giants") were collected [51,52]. The last archaeological trials of these campaigns did not find new evidence and excavations were suspended. Recently, new archaeological and geophysical investigations were carried out from 2013 to 2015 in the framework of a joint project between the Universities of Cagliari and Sassari and the Archaeological Superintendency of Cagliari and Oristano, in which context other significant discoveries were made [53][54][55][56]. Overall, 28 statues representing figures of boxers, archers, warriors, measuring up to 2.20 m in height, 16 models of "nuraghe" and baetyls were pieced together from over 5000 fragments of limestone. Archaeologists identified more than 50 shaft graves. Figure 1a,b indicates the geographic location of the site of study and the area with main archaeological digs carried out in 2015 (Figure 1c). After this campaign, the full extent of the archaeological site appears significantly greater than that supposed by the first researchers in 1970s.

Geological Setting
The geology of the Sinis peninsula, the wide region including the Mont'e Prama site, consists of both sedimentary and volcanic rocks [57]. In particular, the geographic sector closest to prospection area is mainly characterized by sedimentary lithologies. These formations range from messinian fenestrae-rich limestones and breccias of supratidal and intertidal environment with typical benthic fossil content to evaporitic fine-grained limestones, sub-litoral marls and bioclastic limestones (Calcari Laminati del Sinis, Capo S. Marco formation). Quaternary conglomerates, sands and mud deposits, organized in river terraces and alluvial cones, are recognized in the upper part of the sequence.
In the context of 2014 archaeological digs, the stratigraphic sequence in correspondence to the excavations was reconstructed [58]: it is characterized by the presence of geological units related to the sedimentary complex of the Sinis peninsula. In particular, the site presents a thin layer of silt-sandy soils (Pleistocene and Olocene) with a maximum thickness between 40 cm and 50 cm. The soil layer covers a yellowish thick carbonate crust of pedogenetic origin. The hard crust is probably derived from salt precipitation as a result of capillarity action in past climatic conditions and reaches 1 m of maximum thickness. This horizon is observed with good spatial continuity in all over Sinis peninsula and covers previous messinian sediments, mainly composed by carbonate sandstones and marlstones. The area of geophysical surveys is expected to present mean characteristics similar to these results with a certain spatial variability.

ARP Survey
The electrical surveys were carried out using the multi-polar ARP03 device (ARP ® , Automated Resistivity Profiling, Géocarta SA, Paris, France), a measuring system made of four pairs of rolling electrodes arranged in a characteristic V-shaped configuration. Figure 2a shows a picture of the ARP system. The electrodes located on the first pair of wheels represent the current dipole, through which an electrical current of 10 mA is flowing. The other three couples of rolling electrodes are three potential dipoles designed to measure the electrical potential with increasing depth (Figure 2b). The longitudinal distance between the potential dipole electrodes and the current electrodes corresponds to 0.50 m, 1.0 m and 1.70 m, respectively. Moreover, the three potential dipoles constitute of pairs of wheels with axial distances of 0.5 m, 1.0 m and 2.0 m, respectively. Therefore, the resistivity measurements are related to three different depths of investigation, conventionally indicated as 0.5 m, 1.0 m and 1.70 m.
The sensor allows the acquisition of multiple measurements of apparent resistivity through the three pairs of rolling electrodes with a sampling frequency of 22.7 Hz. A radar doppler system is used to measure the distance. Depending on the vehicle speed, raw data can be acquired up to a 10 cm sampling interval.
Moreover, the positioning system, equipped with a differential correction DGPS, enables georeferencing of the acquired data and showing raw resistivity data in real time.
The speed of data acquisition is mainly due to the morphological features of the investigated area (slope, irregular topography, shallow obstacles, plan shape and irregularities, etc.). Usually, the employment of a quad bike to tow the carriage allows a maximum speed of 15 km/h. The acquisition on the ground is typically performed along parallel profiles with variable length. For high-resolution surveys, such as archaeological investigations, the acquisition pattern is composed of profiles positioned at a 1 m distance. The sensor allows the acquisition of multiple measurements of apparent resistivity through the three pairs of rolling electrodes with a sampling frequency of 22.7 Hz. A radar doppler system is used to measure the distance. Depending on the vehicle speed, raw data can be acquired up to a 10 cm sampling interval.
Moreover, the positioning system, equipped with a differential correction DGPS, enables georeferencing of the acquired data and showing raw resistivity data in real time.
The speed of data acquisition is mainly due to the morphological features of the investigated area (slope, irregular topography, shallow obstacles, plan shape and irregularities, etc.). Usually, the employment of a quad bike to tow the carriage allows a maximum speed of 15 km/h. Data were collected in autumn and on a day characterized by optimal weather conditions in order to perform the measurements with humidity conditions of the outermost shallow layers able to guarantee low contact resistances and high-quality measurements.
The acquisition allowed the detection of 277,319 punctual measurements of apparent electrical resistivity along a 27.8 km long route and with a velocity of 7.53 km/h along profiles oriented in an easterly direction. Table 1 summarizes the main acquisition parameters. The spatial distribution of the measurement points is shown in Figure 3a. From this representation, it is possible to observe two distinctive points densities between the two main datasets, partially due to the different acquisition velocities, which are mainly related to the availability of undisturbed wide paths for the towing vehicle but mostly due to the overlap configurations; in fact, the smaller area was surveyed along shorter and orthogonal profiles, even because of superficial archaeological features partially prohibiting vehicle motion, which produced a denser points distribution. While in the northern free area it was possible to perform the acquisition with optimal sampling parameters, for the southern most difficult area, the requirement of spatial coverage and reliability of the local dataset entailed a little higher redundancy of data (rough density of measurements about 29 point/m 2 ).
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 23 The acquisition on the ground is typically performed along parallel profiles with variable length. For high-resolution surveys, such as archaeological investigations, the acquisition pattern is composed of profiles positioned at a 1 m distance.
The ARP survey at the site was carried out on November 21 st , 2014, on a surface of about 2.28 hectares (Figure 1c). Two series of measurements were done: 1. the first dataset was collected in the sector located in the north of the area involved in the archaeological excavations (large 20,155 m 2 ); 2. the second dataset was acquired in the area around the open archaeological digs (2730 m 2 ).
Data were collected in autumn and on a day characterized by optimal weather conditions in order to perform the measurements with humidity conditions of the outermost shallow layers able to guarantee low contact resistances and high-quality measurements.
The acquisition allowed the detection of 277,319 punctual measurements of apparent electrical resistivity along a 27.8 km long route and with a velocity of 7.53 km/h along profiles oriented in an easterly direction. Table 1 summarizes the main acquisition parameters.
The spatial distribution of the measurement points is shown in Figure 3a. From this representation, it is possible to observe two distinctive points densities between the two main datasets, partially due to the different acquisition velocities, which are mainly related to the availability of undisturbed wide paths for the towing vehicle but mostly due to the overlap configurations; in fact, the smaller area was surveyed along shorter and orthogonal profiles, even because of superficial archaeological features partially prohibiting vehicle motion, which produced a denser points distribution. While in the northern free area it was possible to perform the acquisition with optimal sampling parameters, for the southern most difficult area, the requirement of spatial coverage and reliability of the local dataset entailed a little higher redundancy of data (rough density of measurements about 29 point/m 2 ).
The GPS data continuously and automatically acquired during the measurements at 1 Hz sample rate, integrated with isolated points external to the ARP acquisition area, also allowed the creation of a detailed digital terrain model (DTM) of the surveyed area, which provides useful information for the analysis and interpretation of ARP data. Based on DTM results, the site presents a regular topography with a weak slope angle (mean dip angle 4.3°) indicatively oriented from west (top) to east (bottom) for a difference in height of about 16.5 m (Figure 3b). In order to delete both negative resistivity values and outliers, raw data were filtered using a bidimensional median filter operating on a mobile spatial window with a 2 meter radius centered on each measurement point. The filtering algorithm removes resistivity values that are higher than the reference threshold function of the median of values measured within an area centered on every measuring point. Once filtered, the experimental data are interpolated using 2D cubic spline The GPS data continuously and automatically acquired during the measurements at 1 Hz sample rate, integrated with isolated points external to the ARP acquisition area, also allowed the creation of a detailed digital terrain model (DTM) of the surveyed area, which provides useful information for the analysis and interpretation of ARP data. Based on DTM results, the site presents a regular topography with a weak slope angle (mean dip angle 4.3 • ) indicatively oriented from west (top) to east (bottom) for a difference in height of about 16.5 m (Figure 3b).
In order to delete both negative resistivity values and outliers, raw data were filtered using a bi-dimensional median filter operating on a mobile spatial window with a 2 meter radius centered on each measurement point. The filtering algorithm removes resistivity values that are higher than the reference threshold function of the median of values measured within an area centered on every measuring point. Once filtered, the experimental data are interpolated using 2D cubic spline functions over plan surface to permit to obtain a regular grid that can be used to support a successive processing phase or the interpretation of the geophysical anomalies.

GPR Survey
Ground Penetrating Radar methods study both electric and magnetic properties of materials by means of the interaction of microwave signals with the investigated bodies [59]. Electromagnetic waves, transmitted using different superficial probes, can be refracted, reflected and diffracted at the electromagnetic discontinuities and are absorbed passing through the medium. These effects are a function of material properties (e.g., electrical permittivity, magnetic permeability) and signal features (e.g., wavelength, waveform, etc.). These techniques have been widely used in recent decades in engineering, geology and archaeology [60].
In order to validate and interpret ARP data, a comparison with some GPR profiles is proposed. The Ground Penetrating Radar survey was performed in the context of the same research project on a wide area of study, even including the sector investigated by means of the ARP method [31][32][33][34]. The GPR prospection was done using a 15-channels (16 antennas) GPR system (STREAM-X, by IDS), with a 200 MHz central frequency that allows the acquisition of synchronous recordings of 15 parallel radargrams with a fast acquisition speed (about 10-15 km/h). Each radar profile is spaced 0.12 m from the others. Therefore, by means of this instrumental setup, it was possible to collect a big dataset characterized by a high spatial density (horizontal: 12 cm perpendicularly to the moving direction and 9 cm along it: vertical: less than 1 cm). A differential GPS antenna was used for positioning the measurements with a horizontal shift less than 5 cm. This equipment with respect to standard single/dual channels GPR systems includes different advantages such as time efficiency, fixed distance between the fifteen profiles (12 cm), absolute parallelism between all synchronous profiles and 3D highlighting of buried targets.

ARP Results
The geophysical surveys carried out by means of the ARP system permitted to simultaneously obtain three apparent resistivity maps, referring to the three depth intervals. Their measured global resistivity values ranged between 11 and 154 Ω·m, with different responses between each depth map ( Figure 4). First, the shallowest map (Figure 4a), referred to a 0-0.50 m layer, reports the apparent resistivity values ranging from 50 to 117 Ω·m with a data distribution showing maximum frequency at about 80 Ω·m but generally more shifted towards low than high resistivity values for less numerous classes (Figure 4d values ranging from 11 Ω·m to 56 Ω·m with maximum concentration around 25.5 Ω·m and a general symmetrical distribution (Figure 4). The 1.00 m map appears spatially more noisy than 0.50 m one but it is possible to note a good coherence between their resistivity patterns. The deepest map, 0-1.70 m layer (Figure 4c) shows apparent resistivity data ranging from 30 to 154 Ω·m with maximum concentration of around 60 Ω·m, significantly lower than the mean value (92 Ω·m) distributed with a lower spreading (Figure 4f). As is predictable, the results of the 1.70 m map are the smoothest obtained, with clear readability of resistivity patterns and anomalies at a bigger spatial wavelength and a lower presence of both spatially high-frequency signals and noise. Vertical gradients of resistivity ranges are probably related to soil layering properties, in particular, a possible compaction of deepest volumes due to the proximity to bedrock and a moisture gradient for the two shallowest layers due to superficial drying phenomena or different clay concentration with depth. In all the cases, the apparent resistivity patterns show no relation with the topography gradient that characterizes the investigated sector of the site. The result of tillage operations, carried out in the recent past, appear in the form of stripes east-west oriented and in some route traces, are recognizable as double parallel curves on the shallowest map (0.5 m), Figure 5a, b. As is predictable, the results of the 1.70 m map are the smoothest obtained, with clear readability of resistivity patterns and anomalies at a bigger spatial wavelength and a lower presence of both spatially high-frequency signals and noise. Vertical gradients of resistivity ranges are probably related to soil layering properties, in particular, a possible compaction of deepest volumes due to the proximity to bedrock and a moisture gradient for the two shallowest layers due to superficial drying phenomena or different clay concentration with depth. In all the cases, the apparent resistivity patterns show no relation with the topography gradient that characterizes the investigated sector of the site. The result of tillage operations, carried out in the recent past, appear in the form of stripes east-west oriented and in some route traces, are recognizable as double parallel curves on the shallowest map (0.5 m), Figure 5a Looking for archaeological purposes at apparent resistivity patterns, the most significant ones seem to be resistive areas (red tones), which are highlighted in the following maps. Nevertheless, some conductive patterns (blue tones) are also evidenced and reported as signals of potential interest, especially for the top and bottom layers.
In Figure 6, the resistivity map at 0.5 m nominal depth is reported with the indication, bordered in red, of some interesting resistive areas which have vertical continuity even in other two maps. The Looking for archaeological purposes at apparent resistivity patterns, the most significant ones seem to be resistive areas (red tones), which are highlighted in the following maps. Nevertheless, some conductive patterns (blue tones) are also evidenced and reported as signals of potential interest, especially for the top and bottom layers.
In Figure 6, the resistivity map at 0.5 m nominal depth is reported with the indication, bordered in red, of some interesting resistive areas which have vertical continuity even in other two maps. The same areas are reported also at a 1.0 m depth (Figure 7) and 1.7 m (Figure 8). Figure 9 reports the most interesting conductive areas: in particular, the most superficial map (0.5 m, Figure 9a) shows, in its southern part, some conductive regular patterns (indicated as I1), just inside the protected area where digs were concentrated, which could have archaeological relevance such as tombs or similar geometries filled with moist silty or clay soil. The 0.5 m map presents other conductive larger regions which have dimensions not compatible with archeological features and are most probably interpretable as superficial geological variations: this is the case of the G1 region that is confirmed by deeper layers of the dataset and interpreted as a continuous geological body. Indeed, the map at 1.7 m shows wide conductive areas which could have importance to reconstruct paleoenvironment and paleo-landscape of the area (Figure 9b, F3, G3 and H3); area I3 of the same map shows that at higher depths, the conductive anomalies indicated before as I1 inside the protected are less evident and probably their low effect could be due to the fact that the deeper map still includes all the layers between 0.0 m and 1.7 m.

Vertical Analysis of ARP Anomalies and Comparison with GPR Data
To explore the three-dimensional features of the ARP results, in the next paragraph, the data are analyzed by means of their vertical pseudo-sections in correspondence to the most important planimetric patterns. ARP pseudo-sections are also compared with GPR sections acquired over the same paths ( Figure 10). Five vertical profiles have been selected over areas where coverage of both datasets was available: three are located at the north-western part of the ARP-surveyed area, intersecting three resistive anomalies A-B-C in Figures 6-8; the other two profiles are located at the eastern part of the ARP-surveyed area, just outside the protected area partially interested by

Vertical Analysis of ARP Anomalies and Comparison with GPR Data
To explore the three-dimensional features of the ARP results, in the next paragraph, the data are analyzed by means of their vertical pseudo-sections in correspondence to the most important planimetric patterns. ARP pseudo-sections are also compared with GPR sections acquired over the same paths ( Figure 10). Five vertical profiles have been selected over areas where coverage of both datasets was available: three are located at the north-western part of the ARP-surveyed area, intersecting three resistive anomalies A-B-C in  In Figures 11a-c and 12a,b the comparison of radargrams and resistivity pseudo-sections is reported. As also observed in Figure 11a, looking at pseudo-sections, it is possible to notice the difference of apparent resistivity ranges at different depths: mean-resistive data on the top layer; a more conductive layer at middle-depth in which relatively resistive points vertically correspond with the top layer ones; the most resistive layer is the deeper one in which vertical correspondence of relatively high resistivity data with the two top layers is confirmed. This behavior leads to the In Figures 11a-c and 12a,b the comparison of radargrams and resistivity pseudo-sections is reported. As also observed in Figure 11a, looking at pseudo-sections, it is possible to notice the difference of apparent resistivity ranges at different depths: mean-resistive data on the top layer; a more conductive layer at middle-depth in which relatively resistive points vertically correspond with the top layer ones; the most resistive layer is the deeper one in which vertical correspondence of relatively high resistivity data with the two top layers is confirmed. This behavior leads to the consideration that the resistive bodies should have a certain vertical continuity from bottom to top, even in the intermediate layer where a probable higher humidity of soil partially influences apparent resistivity ranges shifting all values towards more conductive ones. The three nominal depths of ARP raw data are indicated with the black horizontal dotted lines on pseudo-sections.
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 23 Figure 11c reports the comparison between the TI060005 radar section and the resistivity pseudo-section at the same position. This line is parallel to the TI040017 radar profile, previously presented ( Figure 11b) and is located about 10-12 m south of this. As Figure 10 shows, this profile intersects the same resistive areas (A and C in Figures 6-8) close to their southern border. The pseudosection shows the same resistive volumes but with slightly lower resistivity values and wider spatial extents. The GPR profile is characterized by the presence of diffuse intermediate signals with two regions where signal amplitude increases: the first one is particularly strong and concentrated between the 60 m and 80 m progressive coordinates, in correspondence with the western resistive body. The second region, located between 125 and 150 m, presents weaker GRP signals than the first one but is also weaker than the corresponding second anomaly in the TI040017 profile. A preferential direction of radar signal is less evident in this case with regards to the TI040017 GPR profile. The depth ranges are the same as in the TI040017 profile for both regions.  the previous one, TI040062. Because of this, anomalous patterns linked to the same body, which are still present, are located at different horizontal coordinates. The deeper ARP anomaly presents the highest resistivity values of all the datasets, saturating the color-scale at 150 Ω·m between 5 m and 45 m, but resistive area can be also considered between 0 m and 55 m. The top resistive area is a bit wider than the range from 20 m to 40 m. Even for this profile, GPR data are in good agreement with the ARP pseudo-section: in fact, a large region with clear sub-horizontal GPR signals is present between 5 m and 40 m, at an approximative depth of 1-2 m.  In Figure 11a, the northern profile (TI040009) is reported. ARP pseudo-section presents two resistive patterns belonging to regions B and C in Figures 6-8. The bigger one, belonging to region B, is located at progressive horizontal coordinates range of 65-90 m and corresponds to parts of the radargram with an intermediate amplitude signal. Moving to the right in the radargram, it is possible to notice the signal weakening in correspondence of conductive areas where microwaves are more absorbed. From 135 m to 140 m, the second resistive anomaly appears as a section of the top linear resistive pattern in region C; at the same coordinates, GPR data present the most noisy traces.
In Figure 11b, the ARP profile and the radargram acquired along the TI040017 line in the south-western side of the study area are represented. The profile crosses two wide resistive anomalies, which are identified as A and C in Figures 6-8. In correspondence of both resistive regions, the radar profile is characterized by the strong amplitudes of the signals. Moving from west to east, a strong reflective region is found from 62 m to 75 m, with a globally complex shape due to the superimposition of multiple diffraction signals. This anomalous region, at depth from 0.5 m to 2.0 m, is spatially related to the first resistive body of the TI040017 pseudo-section, clearly detected also in ARP maps. The second significant radar-reflective region is located from 90 m to 110 m. This signal pattern is related to targets with a decreasing depth from left to right approximately ranging from 1.5 m to 0.5 m. The same trend is observable in the resistivity pseudo-section where the top layer resistive region is right-shifted compared to that of the bottom layer. Figure 11c reports the comparison between the TI060005 radar section and the resistivity pseudo-section at the same position. This line is parallel to the TI040017 radar profile, previously presented ( Figure 11b) and is located about 10-12 m south of this. As Figure 10 shows, this profile intersects the same resistive areas (A and C in Figures 6-8) close to their southern border. The pseudo-section shows the same resistive volumes but with slightly lower resistivity values and wider spatial extents. The GPR profile is characterized by the presence of diffuse intermediate signals with two regions where signal amplitude increases: the first one is particularly strong and concentrated between the 60 m and 80 m progressive coordinates, in correspondence with the western resistive body. The second region, located between 125 and 150 m, presents weaker GRP signals than the first one but is also weaker than the corresponding second anomaly in the TI040017 profile. A preferential direction of radar signal is less evident in this case with regards to the TI040017 GPR profile. The depth ranges are the same as in the TI040017 profile for both regions.
The last comparisons of ARP and GPR data, Figure 12a,b, mainly concern the most resistive and widest area of ARP maps which was indicated with D in Figures 6-8 in the southern part of the prospection. This highly apparent resistivity region is partially located in the restricted archaeological area in the south, and in the bordering part of the free northern area. The two comparative profiles are placed in the free northern area, the first one over the northern limit of the wide resistive pattern and the second one crossing the central higher resistivity part.
In Figure 12a, the first profile, TI040062, is plotted. The pseudo-section is characterized by a large resistive region at its bottom layer, ranging from 15 m to 75 m progressive horizontal coordinate, with the highest apparent resistivity values concentrated between 45 m and 75 m. On the top layer, resistive values are right-shifted and less extended, located between 60 m and 80 m. the GPR profile is in very good agreement with ARP pseudo-section, presenting an intense group of signals between 60 and 80 m with a preferential tilted direction climbing up towards the right. These reflected signals come from variable depths ranging between 0.5 m and 1.8 m. Figure 12b reports the TI100073 comparison profile, which crosses the D in

Discussion
Interpretation of the geophysical results is suggested in Figure 13, where surveyed areas are classified in different color classes by evidence of their archaeological potential: • dark gray color: areas with highly evident ARP anomalies without GPR comparison; • red color: areas with high ARP resistivity contrast to background and consistent high amplitude GPR signals; • orange color: areas with ARP anomalies with good accordance to strong GPR reflections; • yellow color: areas with significant features of ARP anomalies with only partial accordance with GPR signals; • light gray color: areas with relatively weak ARP spatial patterns without clear agreement or comparison with GPR data.
• red color: areas with high ARP resistivity contrast to background and consistent high amplitude GPR signals; • orange color: areas with ARP anomalies with good accordance to strong GPR reflections; • yellow color: areas with significant features of ARP anomalies with only partial accordance with GPR signals; • light gray color: areas with relatively weak ARP spatial patterns without clear agreement or comparison with GPR data. Figure 13. Archaeological interpretation of joint geophysical results in five classes: areas with simultaneous presence of strong ARP and GPR anomalies (R1 and R3a, orange); significant ARP anomalies with locally poor, weak or too noisy, GPR signals (R2 and R3b, yellow); simultaneous presence of strongest ARP and GPR signals (R4a, red); areas with important ARP anomalies without GPR comparison (R4b, dark gray); areas with relatively weak, sometimes coherent, ARP spatial patterns without clear agreement or comparison with GPR data (remaining surveyed area, light gray).
After the joint analysis of ARP and GPR data, it was possible to classify some areas in order of potential archaeological importance. In particular, the most significant areas of the ARP survey could be divided into five classes on the basis of the agreement of apparent resistivity anomalies and the GPR signal behavior or the absence of GPR data. One of the most promising areas is the southern sector internal to the restricted archaeological area, indicated with R4b in Figure 13 (dark gray color, 2730 m 2 ). This complex pattern of conductive and resistive anomalies is characterized by linear and sub-rectangular shapes of dimensions which can be attributed to macro-archaeological features. Nevertheless, the absence of higher resolution data, such as GPR, does not allow to fully analyze tridimensional relationships of the different ARP signals and the sector is classified as very promising but not completely explored with complementary geophysical tools.
The other most important group of anomalies, indicated with R4a in Figure 13 (red color), is contained in a semi-elliptical shape with the main dimensions of 75 m and 25 m (1460 m 2 ), in spatial continuity with the previous anomaly. The area is characterized by a cluster of very high resistivity values with a main dimension of about 40-50 m at the deeper level of ARP measures (150 Ω·m at 1.70 m nominal depth). In this region GPR profiles were able to highlight parts of bigger importance inside Figure 13. Archaeological interpretation of joint geophysical results in five classes: areas with simultaneous presence of strong ARP and GPR anomalies (R1 and R3a, orange); significant ARP anomalies with locally poor, weak or too noisy, GPR signals (R2 and R3b, yellow); simultaneous presence of strongest ARP and GPR signals (R4a, red); areas with important ARP anomalies without GPR comparison (R4b, dark gray); areas with relatively weak, sometimes coherent, ARP spatial patterns without clear agreement or comparison with GPR data (remaining surveyed area, light gray).
After the joint analysis of ARP and GPR data, it was possible to classify some areas in order of potential archaeological importance. In particular, the most significant areas of the ARP survey could be divided into five classes on the basis of the agreement of apparent resistivity anomalies and the GPR signal behavior or the absence of GPR data. One of the most promising areas is the southern sector internal to the restricted archaeological area, indicated with R4b in Figure 13 (dark gray color, 2730 m 2 ). This complex pattern of conductive and resistive anomalies is characterized by linear and sub-rectangular shapes of dimensions which can be attributed to macro-archaeological features. Nevertheless, the absence of higher resolution data, such as GPR, does not allow to fully analyze tridimensional relationships of the different ARP signals and the sector is classified as very promising but not completely explored with complementary geophysical tools.
The other most important group of anomalies, indicated with R4a in Figure 13 (red color), is contained in a semi-elliptical shape with the main dimensions of 75 m and 25 m (1460 m 2 ), in spatial continuity with the previous anomaly. The area is characterized by a cluster of very high resistivity values with a main dimension of about 40-50 m at the deeper level of ARP measures (150 Ω·m at 1.70 m nominal depth). In this region GPR profiles were able to highlight parts of bigger importance inside the wider anomaly which correspond to the most resistive elements (Figure 12a, b). The area is likely to contain archaeological remains, in particular, because of the complex patterns of GPR signals which are difficult to interpret as geological features. The wide extension of this anomaly and its topographic position at the lower altitude of the agricultural field suggest the interpretation as a possible accumulation of stone pieces, even of possible archaeological importance, due to ancient ploughing.
The third most important anomaly is indicated with R1 in Figure 13 (orange color). It is an area 390 m 2 in extent, with 40 m by 13 m main axes, with apparent resistivity values up to 120 Ω·m at 1.7 m nominal depth. After the integrated analysis of ARP and GPR data (Figure 11b,c), the anomaly could be interpreted, even for this case, as a localized concentration of stone bodies and fragments. In fact, GPR reflections and diffractions are characterized by high amplitude and spatial concentration without lateral continuity, at least at depths with an acceptable Signal to Noise ratio. Similar properties are recognizable for area R3a, which has an extension of 380 m 2 , 16 m wide and 38 m long. The comparison of ARP and GPR data (right anomaly in Figure 11b,c) emphasized a correspondence, even if weak, between the two datasets. Resistivity values are on the same order as anomaly R1 (up to 120 Ω·m) and GPR signals present high amplitude in the northern section of the anomaly but are still readable in the southern one. Especially at the northern section, GPR anomalies do not show lateral continuation outside the resistive volumes. Because of these considerations, anomalies in areas R1 and R3a can be potentially considered for fruitful archaeological investigations.
To interpret the other two regions bordered in Figure 13, R2 (187 m 2 , 22 m × 12 m) and R3b (284 m 2 , 47 m × 8 m), yellow color, the joint analysis of ARP and GPR data was less effective, because in correspondence of the area R2, radar signal was present but not intense as in the other cases, while in correspondence of norther section of area R3b, it was affected by the high noise level. Nevertheless, lateral apparent resistivity contrasts and spatial features suggest the opportunity for further investigations with other complementary geophysical methods and archaeological trials.
The planimetric analysis of apparent resistivity data in the remaining areas (light gray color in Figure 13, 17,830 m 2 ) evidences a complex scattered distribution of spatial patterns with a low resistivity contrast which could be indicative of the presence of in situ archaeological features, less sensitive than possible accumulation of stone fragments to volume analysis, such as ARP survey.
Aiming at the estimation of the effective investigation depth of the ARP survey, an inversion procedure, based on a tetrahedral Finite Elements structure, was performed on two apparent resistivity profiles by means of the commercial software ERTLab. Pseudo-sections TI040017 and TI040062 were selected by considering the good concordance between apparent resistivity anomalies and GPR signals. The resistivity models were set up with 3D cells of 20 cm edge dimensions and inverted with smoothness constraints to achieve the final resistivity models. The tomographic results confirm the location of resistive bodies and the presence of an intermediate conductive layer (Figure 14). Some artifacts are present due to the scarcity (only three levels) of data distributed across the pseudo-sections: this is most evident for a shallow layer over the top of the first apparent resistivity level where most resistive values are distributed and characterized by a very high discontinuity along the horizontal profile. The same effect is quite common, even for traditional acquisition ERTs, typically with minor thicknesses [61]. The deepest bodies under the intermediate conductive layer still have a resistive behavior, but values are significantly lower than in the case of pseudo-sections. the wider anomaly which correspond to the most resistive elements (Figure 12a, b). The area is likely to contain archaeological remains, in particular, because of the complex patterns of GPR signals which are difficult to interpret as geological features. The wide extension of this anomaly and its topographic position at the lower altitude of the agricultural field suggest the interpretation as a possible accumulation of stone pieces, even of possible archaeological importance, due to ancient ploughing.
The third most important anomaly is indicated with R1a in Figure 13 (orange color). It is an area 390 m 2 in extent, with 40 m by 13 m main axes, with apparent resistivity values up to 120 Ω·m at 1.7 m nominal depth. After the integrated analysis of ARP and GPR data (Figure 10b, c), the anomaly could be interpreted, even for this case, as a localized concentration of stone bodies and fragments. In fact, GPR reflections and diffractions are characterized by high amplitude and spatial concentration without lateral continuity, at least at depths with an acceptable Signal to Noise ratio. Similar properties are recognizable for area R3a, which has an extension of 380 m 2 , 16 m wide and 38 m long. The comparison of ARP and GPR data (right anomaly in Figure 10b, c) emphasized a correspondence, even if weak, between the two datasets. Resistivity values are on the same order as anomaly R1a (up to 120 Ω·m) and GPR signals present high amplitude in the northern section of the anomaly but are still readable in the southern one. Especially at the northern section, GPR anomalies do not show lateral continuation outside the resistive volumes. Because of these considerations, anomalies in areas R1a and R3a can be potentially considered for fruitful archaeological investigations.
To interpret the other two regions bordered in Figure 13, R2 (187 m 2 , 22 m × 12 m) and R3b (284 m 2 , 47 m × 8 m), yellow color, the joint analysis of ARP and GPR data was less effective, because in correspondence of the area R2, radar signal was present but not intense as in the other cases, while in correspondence of norther section of area R3b, it was affected by the high noise level. Nevertheless, lateral apparent resistivity contrasts and spatial features suggest the opportunity for further investigations with other complementary geophysical methods and archaeological trials.
The planimetric analysis of apparent resistivity data in the remaining areas (light gray color in Figure 13, 17,830 m 2 ) evidences a complex scattered distribution of spatial patterns with a low resistivity contrast which could be indicative of the presence of in situ archaeological features, less sensitive than possible accumulation of stone fragments to volume analysis, such as ARP survey.  The spatial distribution of inverted data substantially illustrates the same patterns of pseudo-sections but most resistive values were obtained on the top soils instead of the bottom ones where the most intense resistive anomalies were located in the pseudo-sections: this opposing behavior was probably caused by data density and distribution which bring ERTs to very smooth variations at middle-bottom depths and to apparently noisy variations on the top depths. Comparing the ERT results and pseudo-sections, it is possible to observe a good agreement of depths for apparent interfaces and spatial resistivity transients, with a slight overestimation of depths under conductive volumes. A similar behavior is also described by Papadopoulos et al. [39,40]. Combining ARP, 0.50 m depth, with DTM data (Figure 15), it is possible to evaluate the independence and spatial relations of the widest ARP anomalies to each other. No clear spatial link is evidenced between the most promising areas identified in Figure 13. The only weak relationship which could be notice is on the eastern bottom areas where the strong linear anomalies are in connection with a northern wide area characterized by less intense resistive anomalies (D1 and E1 in Figure 6). The same connection is no longer present at deeper depths where the northern area (E2 in Figure 7 and E3 in Figure 8) gradually becomes less resistive and disconnects from the strongest anomalies to the south. This feature, although the low contrast of apparent resistivities, could be a marker of a geological transition in the investigated volumes, rather than archaeological remains. From the same three-dimensional comparison, no connection is recognizable along the eastward direction between resistive bodies under areas located at different altitudes. Aiming at the estimation of the effective investigation depth of the ARP survey, an inversion procedure, based on a tetrahedral Finite Elements structure, was performed on two apparent resistivity profiles by means of the commercial software ERTLab. Pseudo-sections TI040017 and TI040062 were selected by considering the good concordance between apparent resistivity anomalies and GPR signals. The resistivity models were set up with 3D cells of 20 cm edge dimensions and inverted with smoothness constraints to achieve the final resistivity models. The tomographic results confirm the location of resistive bodies and the presence of an intermediate conductive layer (Figure  14). Some artifacts are present due to the scarcity (only three levels) of data distributed across the pseudo-sections: this is most evident for a shallow layer over the top of the first apparent resistivity level where most resistive values are distributed and characterized by a very high discontinuity along

Conclusions
A geophysical prospection of a prehistoric archaeological site in Sardinia, Italy, is proposed. In particular, results from an Automated Resistivity Profiling (ARP) survey over a 2.28 hectares area are analyzed and discussed in comparison with Ground Penetrating Radar (GPR) sections over the same area. The purpose of the prospection was to reveal the archaeological potential of a wide area surrounding a site accidentally discovered and partially excavated in the 1970s when a necropolis and related archaeological findings (e.g., fragments of large stone statues and funeral monuments) were found. The site was abandoned until the joint archaeological-geophysical scientific project (2013)(2014)(2015) in which the ARP and GPR prospections were carried out.
The ARP survey covered a wide area localized at the northern side of the historic archaeological site. Apparent resistivity data were simultaneously acquired at three different depths of investigation (0.5 m, 1 m and 1.7 m) and exhibit a wide range of variability (from 15 Ω·m to 250 Ω·m). The most intense planimetric patterns of the apparent resistivity maps were characterized by wide resistive regions with irregular shapes and linear spatial extensions of the order of tens of meters, but even some conductive anomalies were found which were characterized by sub-regular shape and dimensions of possible archaeological interest. Some small apparent resistivity patterns presented spatial coherence but weak local resistivity contrast; therefore, it is difficult to assign them to archaeological features without more detailed measurements of the possible artifacts.
In order to verify preliminary interpretations, the tridimensional distribution of apparent resistivity measures was compared with high resolution GPR data acquired along linear profiles. In most cases, the results show a good correspondence between apparent resistivity anomalies and GPR signals. Based on the agreement between ARP and GPR datasets, it was possible to define the surveyed area in five classes which varies in function of the two kind of signal intensities and their complementary spatial features. Two very promising and neighboring wide areas were delimited for an approximate total extension of 4190 m 2 in the southern part of the surveyed region, corresponding to the limit of the archaeological restricted area (inside and outside). GPR comparison was possible only for one of them (and was positive), but their continuity and ARP data parameters (intensity and spatial features) are evidence of their global archaeological potentiality. Of the other four large resistive regions (in the range of 190-390 m 2 ), two were characterized by the spatial correspondence of strong ARP and GPR anomalies, while the others, presenting high apparent resistivity contrast and shapes potentially attributable to non-geological origin, presented GPR signal of medium amplitude or locally covered by strong noise. Consequently, they were classified as susceptible in the first case to include buried archaeological evidence, and in the second case, suggested further investigations with other complementary geophysical methods and archaeological trials. The remaining wide area (17,830 m 2 ) showed a complex scattered distribution of planimetric patterns with low resistivity contrast which could be indicative of the presence of in situ archaeological features, less sensitive than possible accumulation of stone fragments to volume analyses, such as ARP.
In order to validate ARP apparent resistivity patterns and nominal map depth, two sample pseudo-sections were inverted, substantially confirming their validity but obtaining few differences in terms of hierarchy of most resistive values between top and bottom subsoil and partially in terms of depths in correspondence with the most conductive layers. The depth values of the pseudo-sections and the ERTs bottom interfaces are still compatible with the geological investigations summarized in Section 2.2 considering the ordinary spatial variability of geological features.
In the present case study, the ARP method was confirmed as a reliable and expeditious geophysical technique able to investigate wide regular surfaces and to identify areas likely to contain archaeological features. Thanks to the simultaneous acquisition of positioning data by means of a differential GPS system mounted in the middle of the electrical array, ARP data were instantly geopositioned and an accurate georeferencing was automatically performed. Very good performances were evidenced in indicating large subsoil volumes with a high apparent resistivity contrast. Some spatially coherent patterns were also recognizable with weak apparent resistivity contrast for which further investigations could be useful. Overall, ARP proved to be an effective geophysical technique for the preliminary assessment of the archaeological potential of investigated areas and to direct more concentrated, expensive and time-consuming higher resolution prospections. Furthermore, joint interpretation of ARP and GPR measurements could be used to differentiate potential archaeological zones, demonstrating the effectiveness of the integration of both extensive non-destructive methods to address preliminary researches concerning archaeological studies.
Funding: This research was funded by Regione Autonoma della Sardegna, LR 7/2007. The APC was funded by Fondazione di Sardegna.