Exploring the Consistency of Data Collected in Archaeological Geophysics: A Case Study from the Iron Age Hillfort of Villasviejas del Tamuja (Extremadura, Spain)

: Different geophysical methods applied at the settlement of Villasviejas del Tamuja (Botija, Spain) have identified robust anomalies located at the same position, but some anomalies are reflected by only one method. Furthermore, analysing the spatial correlation of these anomalies is of fundamental importance for obtaining a correct archaeological interpretation. In this work, we analysed the main results of electrical resistivity tomography (ERT), ground-penetrating radar (GPR) and magnetic gradiometry methods in a particular area of the archaeological site. In this analysis, we performed graphical and numerical spatial correlation analyses of the anomalies and observed strong agreement among the results provided by each method. Certain anomalies were reflected only in the magnetic and ERT studies. The results highlight the importance of applying several geophysical methods and performing spatial correlational analyses. Furthermore, the methodology that we have applied to evaluate the spatial correlation offers interesting results. R.J.O., C.C. and J.F.B.; software: C.P., B.C., C.C. and R.J.O.; validation, C.P. and B.C.; formal analysis, C.P. and B.C.; investigation. C.P., B.C., M.T.d.T., C.C. and R.J.O.; resources, C.P., B.C., C.C., R.J.O., and V.M.; data curation, C.P., B.C., C.C. and R.J.O.; writing—original draft preparation, C.P.; and editing, C.P., B.C., R.J.O., J.F.B. and visualization,


Introduction
The archaeological site of Villasviejas del Tamuja is located in the municipality of Botija (Province of Cáceres, Spain); it is a fortified settlement that was founded at the beginning of the 4th century BC, and it remained active until the 1st century BC. Two walled enclosures (A and B in Figure  1) reinforced by artificial ditches and natural defences provided by the Tamuja River and the Verraco stream encompass a total area of approximately 7 ha.
From a geological point of view, the study area is part of the domain of vertical folds of the Central Iberian Zone of the Iberian Massif [1], where the upper Pre-Cambrian successions occupy extensive outcrops. These are metasedimentary series composed of schists and metasandstones, called the Schist-Greywacke Complex (SGC) [2] (Figure 1). The SGC has been divided into different lithostratigraphic units corresponding to the area of study by the Domo Extremeño Group [3]. The sediments are interpreted as having been laid down in submarine fans, slopes and channels [4,5]. The series are characterised by low-grade regional metamorphisms and are affected by different The settlement of Villasviejas sits on a promontory consisting of metasedimentary rocks corresponding to the SGC unit, mainly slates. A geological study carried out at the settlement has allowed the geological characteristics both on and below the surface to be determined [8]. The slates are affected by primary schistosity with differentiated highlights. This structural disposition influenced the initial construction of the settlement and the development of the defensive wall. The slate substrate is located at a depth between 1 and 3 m, and beginning at an approximate depth of 7 m, these slates are affected by a possible granite intrusion that would have generated a bulge and a contact aureole leading to compositional disagreement with the host slate. The effects on the slates have been observed in the potholes [9] in the bed of the Tamuja River that delimits the hillfort to the west.
Villasviejas is an important archaeological site for the study of the Second Iron Age and the romanization process in central Iberia, because it is one of the enclaves of this type of greater extension and monumentality within the region, and it is one of the best-known study cases. Therefore, it serves as a reference for knowledge of the emergence of the hillforts, the first great human settlements in the area before the spread of Roman towns.
This type of settlement reveals the great demographic, economic and social changes of the mentioned period, considering them as the most immediate precedent of the urban phenomenon in this part of Europe. However, knowing the anatomy of these sites-their internal organization, urban planning, domestic architecture and defensive structures-represents a challenge if we depend only on conventional excavation work. In addition to the high cost in time and resources, large excavations pose conservation problems that make the management of these large and complex archaeological areas unsustainable. Bearing this in mind, we approach the study of these sites through the integrated use of a series of prospecting and detection methods that allow obtaining abundant information while minimizing destructive tasks.
To understand this spatial organization and diachronic evolution of the hillfort, we carried out a detailed multidisciplinary study at the settlement that covered the application of geophysical methods, geological and archaeological characterization and digital cartographic representation [10]. We elaborated LiDAR digital terrain model (DTM) hillshade and simple local relief model (SLRM) derived from the LiDAR DTM ( Figure 2). SLRM [11,12] is a trend model useful to enhance the visibility of small-scale features; Figure 2 shows blue surfaces with a relatively lower position than the ones in red, that show a relatively higher position. Unfortunately, this visualization technique does not detect any of the structures identified with geophysics. We also used the information provided by previous excavation efforts ( Figure 2) that have been carried out since the mid-19th century by a different research team [13] for comparison of the results. The aim of this paper is to present the results obtained from a detailed study in a particular area of the archaeological site ( Figure 2) by applying electrical resistivity tomography (ERT), groundpenetrating radar (GPR) and magnetic gradiometry methods. The archaeological excavation is still pending. The study zone covers an area of 20 × 9 m and offers optimal conditions for geophysical work, since the topography is smooth, and the terrain is free of brush and obstacles such as trees or large stones. We detailed the geological knowledge of the substrate in this zone to support the interpretation of geophysical methods. Currently, it is very common in archaeological studies to apply several geophysical methods to investigate anomalies that correspond to the same structure, because different methods may reveal different aspects of the subsurface [14][15][16][17]. In this study, we analysed the main anomalies revealed by each method and evaluated the graphical and numerical spatial correlation.
The application of geophysical methods for mapping cultural remains first emerged in the United States and United Kingdom between 1920 and 1930, with the use of electrical and magnetic methods [18]. These methods have since been widely applied because they allow knowledge of the subsurface to be obtained more quickly than is possible through archaeological excavation, and thus are of great use for defining the study location and planning all further actions to be performed [19]. Currently, there is a great diversity of geophysical methods that can be applied to an archaeological site, such as GPR, vertical magnetic gradiometry and ERT.
With the emergence of these techniques, it has become possible to gain in-depth knowledge of a location; however, there has not always been an in-depth study of the acquired data, as these data frequently show an apparent lack of information. In many cases, this occurs due to a low signal-tonoise ratio. Noise in large quantities masks useful information, resulting in low perceptibility of buried structures. Currently, in-depth study using advanced mathematical techniques is a research area that is undergoing extensive development because of its high applicability in important contexts such as medical and digital image processing. Even data of very low quality can be used as input data to produce images of better quality. Image fusion algorithms are one example of this type of processing [20].

Magnetic Gradiometry
In magnetic methods, variations in the intensity of the Earth's magnetic field are measured to probe the different degrees of magnetization of the materials found in the subsoil. These differences in the magnetic field depend on both induced and remnant magnetizations. These methods enable the exploration of large surfaces within a relatively short time; but in general, without knowledge of the geometry, the depth of structures cannot be established with precision. Therefore, such a method is suitable to carry out a general exploration of a site to preliminarily identify the zones of greatest interest.
In this study, we used a Grad601-2 fluxgate gradiometer from Bartington [21,22] in grid mode with an interval of 50 cm between the data lines and 8 observations per meter. The data were initially processed using Terra Surveyor software to correct errors derived from acquisition, such as stripping (imbalanced between sensors) and staggering (mismatch in horizontal displacements), to perform filtering for the correction of noise problems and other operations such as field transform derivatives to enhance possible alignments. Once the data was exported to a geographic information system (GIS) environment and vector digitized into a polygon geometry, we explored the data distributions and maximized the detection of traces of archaeological interest.
With the aim of analysing the anomalies revealed by the magnetic method in detail, we carried out a pseudo three-dimensional (3D) study based on a grid composed of 19 parallel profiles of 20 m in length and spaced at 0.5 m, with an interval between electrodes of 0.5 m, which yielded results with a satisfactory resolution. The whole area studied was 20 × 9 m, and the greatest depth reached in the study was 4 m. We used an ABEM Terrameter LS resistivity meter with 42 electrodes and multi gradient array. With this configuration, different sets of measurements were performed with the current and potential electrodes at different locations, enabling very stable field data acquisition with a good signal-to-noise ratio [27]. The starting point for each profile was located in the south and the locations of the sampling points were established by using a differential Global Positioning System (GPS) receiver (Trimble R8) with Real Time Kinematic (RTK) corrections received from the Extremadura Positioning Network. The data were georeferenced and processed using the RES2DINV and RES3DINV inversion programs from GeoTomo Software based on the least-squares method involving the finite-element and finite-difference methods.

Ground-Penetrating Radar (GPR)
A GPR device operates on the basis of electromagnetic (EM) pulse waves transmitted by an antenna (Tx) into the subsurface. When radar pulses encounter materials with different dielectric properties, part of their energy is reflected by these materials and returns to the surface where it is collected by the receiver antenna (Rx) if certain theoretical physical constraints are respected [28,29]. The simple model that explains the physical conditions that must hold to enable object detection using this method is described by the Fresnel zone equation [30] where ∆l is the minimum size of an object that can be detected at a distance r from a GPR antenna emitting electromagnetic pulses with a central frequency fc that propagate through the soil with a velocity v. In a geological material = √ ⁄ , where is the dielectric constant. The resulting data are represented in the form of radar sections (radargrams), in which the xaxis represents the antenna survey line and the y-axis represents the propagation time of the radar wave on a given path. The amplitude at each time is related to the electromagnetic properties of the structures by which the waves are reflected (electrical resistivity, dielectric constant = relative permittivity and magnetic permeability). The clearest reflections are produced when the radar pulses are reflected at the interface of two materials with significantly different dielectric constant values. Smaller amplitudes are produced by reflections at the interfaces of materials with dielectric permittivities that are only slightly different. The interpretation of the reflection profiles enables the estimation of the locations of buried structures and their dimensions in the vertical plane. Many radar profiles collected from adjacent locations can be merged to create a three-dimensional map of reflection amplitudes, which may be displayed using 3D visualization techniques either as isosurface renderings or as 2D slices (cutting planes in 3D) [31].
In the GPR survey carried out at Villasviejas del Tamuja, a Geophysical Survey Systems, Inc. (GSSI) equipment was used, including a SIR-3000 GPR controller with a 400 MHz monostatic antenna mounted on a cart with a calibrated tachometer. The survey area of 25 × 34 m was covered by 51 parallel transects with a spacing of 0.5 m, and data were collected in a zig-zag fashion. Table 1 presents the configuration used for the GPR measurements.

Magnetic Gradiometry
The anomalies were classified into a series of discrete categories to identify wall foundations (certain or probable) and thermoaltered surfaces ( Figure 3). Within the latter group, we differentiated among possible pavements, burnt structures and possible combustion structures (kilns and fireplaces).
Based on the geological characteristics of the environment, the soil appears to present mediumlow magnetic properties, which are, however, enhanced by the presence of basic igneous rocks and a likely high concentration of iron oxides in the surface horizon derived from anthropogenic activities. The result is a soil with positive magnetic values, as shown by the magnetic readings taken in the field. Therefore, the results of the geomagnetic survey depend on the contrast between this soil and the construction materials used.
Negative gradient values correspond to walls built with rocks in the area. These rocks are likely to be igneous and metamorphic with relatively low to medium susceptibilities; they are of diamagnetic and induced origin. They correspond to the local materials used in the construction of the wall basements that are already known from the excavated areas of the site, such as slates and granites. The identified walls were isolated in the interpretation as a separate category ( Figure 3D).
Positive gradient values correspond mostly to fillings with materials that exhibit medium magnetic susceptibility, such as the remains of materials subjected to prolonged medium-high temperatures, fragments of pottery or construction materials such as bricks. In some cases, the remains of compacted earthen pavements containing organic material, ceramic fragments and clay subjected to heat may be present.
Dipolar anomalies correspond to two small adjacent areas, one with positive and the other with negative readings of diverse origins; some indicate the presence of ferromagnetic elements, such as metal remains, resulting from anthropogenic activities. When possible, focal dipolar anomalies were distinguished, which correspond to the presence of isolated metals or focal deposits that could have an archaeological origin. Other dipolar anomalies of high intensity and undefined geometry were interpreted as the result of combustion. Finally, a smaller number of dipolar anomalies were classified as indeterminate, as we were unable to specify their origin.

Electrical Resistivity Tomography (ERT)
First, we processed each profile separately to identify the main anomalies and study their continuity along the different profiles; the inversion results had an error of less than 2%. The results are shown in Figure 4, in which three main anomalies that are present in almost all the profiles can be seen. The first anomaly (A1) has a resistivity of 150 Ωm, is located close to the surface and is approximately 0.5 m wide. Considering that we performed the work in February following rainy days, A1 corresponds to filler materials (sediments and archaeological deposits) with high water contents. The second anomaly (A2), with a resistivity equal to 700 Ωm, is present in profiles P1, P2, P3, P6 and P7 at a horizontal distance of 4-7 m from the origin and at a depth between 1 and 2 m. This anomaly is also present in P17 and P18 at 6-10 m from the origin and in P19 at 8-12 m. The third anomaly (A3) has a resistivity of 900-1000 Ωm and appears in P4, P5, P14, P15 and P16. The position and dimensions of this anomaly are similar to those of A2; hence these anomalies could belong to the same structure. Another anomaly with lower resistivity values than the surrounding medium, called A4, is present in the P5-P13 profiles, more clearly evident in the results from the 3D study. We also carried out 3D inversion and obtained the distribution of the resistivity in a cube with dimensions of 20 × 9 × 4 m. From this distribution, we obtained slices at different depths ( Figure 5) that reflect the continuity of the anomalies. At depths of 0.4 and 0.7 m, and with less significance at 1.5 m, we observed the highest resistivity values along lines with NS and EW orientations. Considering the results obtained by applying the ERT method in several excavated zones at this archaeological site [10], we concluded that these features correspond to walls. These resistivity values correspond to the A2 and A3 anomalies. Furthermore, another anomaly appears with a resistivity value of 250 Ωm and a circular shape from 1.5 to 3.8 m in depth; this anomaly is located close to the centre of the grid. This anomaly could correspond to that labelled as A4 from the P5-P13 profiles in Figure 4, although we observe from those profiles that the depth reached by this anomaly is shallower than that indicated by the slice representation; additionally, the anomaly is not so clearly defined. Nevertheless, the anomaly also appears in the magnetic study ( Figure 3). To analyse these anomalies in detail, we generated 3D figures with the corresponding resistivity values. To generate Figure 6a, we selected the A2 and A3 anomalies, and the results show the walls discussed above. Additionally, to generate Figure 6b, the A4 anomaly and the closest and surrounding resistivity values were selected, and the results show the structure that corresponds to the circular anomaly, presented in Figure 5. The two different structures are shown together in Figure  6c.

Ground Penetrating Radar (GPR)
GPR data processing was performed using the commercial software RADAN-7 (GGSI, Inc.) in two phases. First, some transects were individually processed to define the processing stages and their parameters were set to produce well-defined reflection profiles. As the quality of the signals acquired in this study was good, conservative flow processing was applied (Table 2), mainly with the aims of increasing the signal amplitude (gain correction) and, more importantly, increasing the vertical and lateral resolutions (deconvolution operation). The processing pathway defined in the parametrization phase was then applied to all transects, and a 3D amplitude map was obtained by averaging the amplitudes of the reflection profiles [32].   Figure 7 shows a map on which all GPR profiles are marked and the subset of them used in this work is highlighted by the smaller rectangle. The portions of the 18 radargrams that correspond to the analysed area are also presented in Figure 7. From initial observation of this set of profiles, it is evident that many well-defined reflections are present in almost all of them, at different scales of intensity and size. The majority of these reflections are located in the soil layer between the surface and a depth of 1.5 m. These features are characteristic of the typical GPR reflection scenario observed at well contrasted archaeological sites. A closer analysis of the sections allows us to conclude that some of these reflections continue into adjacent profiles. The common marks across adjacent radargrams reveal some extended reflectors crossed by multiple sections and provide an idea of the lateral distribution of the buried structures. These observations are compatible with alignments produced by building remains.  Figure 8 (left) shows a horizontal slice at 0.7 m ± 0.1 m in depth taken from the 3D amplitude reflection map computed from the GPR dataset after processing. As was already suspected from the observation of the radargrams presented in Figure 7, this slice clearly defines a set of alignments compatible with building remains. The area under analysis in this work is marked by the blue rectangle, and four of the transects that define the image are marked with green lines on the slice and are shown in the right panel. On each profile the positions of some high-amplitude reflections that are identified on the related transects are indicated. We note that the perturbations in the travel time caused by topographic variations are very small (~1°) and do not need correction. The 3D amplitude map computed with RADAN-7 was exported in ASCII format using a rectangular grid with mesh dimensions of dx = 0.2 m, dy = 0.06 m and dz = 0.1 m. This process was accomplished by resampling the image to obtain a cubic mesh of square cells with edges of 0.15 m and applying 3D median filtering to obtain similar spatial frequency content.
A two-fold approach for visualizing the 3D GPR results was used. First, a 3D image of the reflecting objects was obtained by rendering a selected amplitude range and assigning transparency to the remaining amplitudes ( Figure 9). In the second step, horizontal slices of the absolute amplitude, corresponding to stacking of the energy at a depth of Z ± 0.1 m, were represented using the same grid cell sizes (dx and dy). Figure 10 displays six of these horizontal cuts from depths that coincide with the ERT model.

Anomaly Correlation
The relative effectiveness of the three geophysical methods used in this case can be determined in a decisive way based on the characteristics of the features that are the target of the exploration. This approach refers to the archaeological stratification (thicknesses and types of sediments) and to the depth and composition of the structures. Regarding the former, the investigated area shows increasing sedimentation in the S-N direction due to the steep slope towards the boundary of the walled enclosure. However, the volume of sediments is generally small and decreases towards the southern part of the area. The fill that covers the spaces between the structures has a clear signal since the building material of the walls is mostly mud bricks and rammed earth. In general, the foundations of the structures were built by means of stone masonry using local materials (slates). Since the foundations appear at very shallow depths (sometimes only approximately 20 cm from the surface), they show the maximal contrast with the surrounding materials. For this reason, although the area has not been excavated, we can obtain information about the structures that are present.
As noted previously, we observed two different types of anomalies by applying the ERT method, one corresponding to walls, and the other corresponding to a circular structure. As the magnetic and GPR methods also both reveal clear anomalies corresponding to walls, we performed a correlation analysis only for this anomaly type. In particular, we conducted the analysis based on the results corresponding to a depth of 0.7 m because at this depth, the anomalies are best revealed. Figure 11 shows the anomaly distribution obtained from the magnetic method (Figure 11a), ERT ( Figure 11b) and GPR (Figure 11c) at a depth of 0.7 m. We observe that although the anomalies corresponding to walls are clearly detected by all three methods, the locations and dimensions are not shown as clearly by the magnetic method as by the other methods. In general, we observe that the zones that are devoid of anomalies coincide among the three methods, but some elements that are evident in the results of one method do not appear in the others. For that reason, to conduct a detailed analysis of the results, we carried out a spatial correlation analysis. First, we isolated the anomalies corresponding to walls ( Figure 12) from each set of data using the same grid spacing and the same number of nodes (a total of 5865). Then, we developed a binary system by assigning a value of 1 to the nodes where anomalies were observed and a value of 0 to the rest to allow a numerical correlation analysis to be carried out by multiplying these values. In this way, at nodes where an anomaly is present, we obtained a value of 1, and otherwise, we obtained a value of 0. We also calculated the corresponding values of the Pearson correlation coefficient (Table 3), which measures the statistical relationship between two variables. In particular, this coefficient is appropriate for geophysical studies applied in archaeology to identify any correlations existing between datasets collected using different geophysical methods [15,33,34]. The obtained values of this coefficient show that the best correlation is between the ERT and GPR results and that the worst correlation corresponds to the GPR and magnetic results. To analyse the spatial distribution of the correlations between the anomalies, we show the nodes at which the identified anomalies appear for each method in Figure 13. In Figure 13a, the positions at which the anomalies obtained via ERT are located are represented in yellow, the positions at which anomalies are obtained via GPR are represented in blue, and the positions at which the two types of anomalies coincide are represented in grey. Moreover, the magnetic anomalies are shown in red, and the nodes at which there is concordance between the ERT and GPR anomalies are shown in orange and purple, respectively. Finally, the points at which all three anomaly types coincide are shown in brown. Figure 13b-e shows the correlations for all possible pairs and for all three types of anomalies together. These figures reveal that the areas that are devoid of anomalies are consistent among all three methods across most of the study area. We also calculated the percentages of the nodes for which the different types of anomalies coincided (Table 4). These results are fundamentally the same as those shown in Figure 13, with the best correlation being between the ERT and GPR approaches; however, the specific values allow us to quantitatively compare the areas of coincidence.

Discussion
Although the archaeological exploration of Villasviejas dates back to the mid-19th century, the systematic research that began in the 1970s allowed discovery of the remains of many building structures. The excavation findings served as useful references to test the validity of the survey results (the area studied in this work has not yet been excavated). A clear coincidence was observed regarding the shape and orientation. Additionally, results from magnetic survey next to excavated areas allowed us to roughly assign the ranges of geomagnetic survey values corresponding to these structures. These clear results, in combination with the high contrast of the ERT and GPR data, made it possible to go beyond a mere quantitative correlation between data from different methods and provided some clues regarding the structures identified.
It is clear that in general, the three geophysical methods used in the study area offer a fairly good comparison for the identification of archaeological targets. Notably, as has been shown by the results of the quoted excavations in other areas at the archaeological site, most of the wall foundations are very shallow, sometimes almost at the surface. Hence, there is no attenuation of the sensor response.
Additionally, the building technique used for the structures, which consist of drystone foundations as a basement for walls made of earth, helps to facilitate anomaly detection. When these structures collapse, the filling deposits within the rooms are very clean stone blocks, resulting in a very "clean" and homogeneous sensor response.
From the graphical and numerical correlations among the three types of anomalies, we observe that the spatial concordance of the three methods is poor: the greatest coincidence is found for the GPR and ERT results, with the greatest area of positive correlation. This result is also supported by the Pearson coefficients, for which the highest value corresponds to the ERT and GPR results (Table  3). Similar results were obtained by Kvamme [15] on the basis of the local Pearson coefficients in small neighbourhoods; however, this author obtained a better correlation for the magnetic and GPR results than for the magnetic and ERT anomalies. Meanwhile, in this work, we found that the magnetic results best match the ERT results (Table 3 and Figure 13).
Moreover, Figure 13a indicates that the areas that are devoid of anomalies are in agreement among the different methods, except in the northern part of the study area, where only the magnetic method seems to reflect a structure. The same results are shown in detail in Figure 13b-e, where we also observe a correlation for zones without robust anomalies.
Therefore, the visual impression, which is confirmed by statistical validation, suggests that while the magnetic method can record subtle variations in the composition of buried deposits, ERT and GPR can reveal high-definition anomalies corresponding to structures.
The cylindrical anomaly revealed by the ERT data is not seen in the GPR data, probably because the corresponding resistivity is lower than the resistivity of the walls; only a small discrepancy with the surrounding materials is observed, and this difference is difficult to detect with certain methods. Considering the topographic setting of the building, which is on a steep slope facing the walled enclosure of the settlement, it is possible that the hillfort could have a terraced surface with a level below the ground. The artificial division of the rooms in the deep ERT readings supports this idea. In this context, this anomaly could correspond to a cylindrical structure designed for storage or a pit excavated in the bedrock. Bearing in mind that the ERT study was carried out in February after rainy days, this structure would be occupied by sediments with a high water-content and therefore display a low electrical resistivity. The magnetic data also show this structure, which appears as a positive anomaly, corresponding to the filling materials.

Conclusions
We analysed the results obtained by applying magnetic, GPR and ERT survey methods over an area with dimensions of 20 × 9 m. These three methods reveal the main anomalies, but we observe that some structures are not reflected by all of these methods. This finding reveals the importance of applying different geophysical methods to enable a complete interpretation of anomalies. With magnetic gradiometry, it is possible to cover a wide area within a short time, thereby obtaining the general locations of the main anomalies with variations in details based on depth. The ERT method can yield detailed images (2D and 3D) of the structures, although the field work for this method takes a long time. Finally, the GPR method is quick to implement in the field and provides 2D and 3D underground images, but it does not reflect structures with low electrical resistivity values. Regarding the duration of the acquisition work, in this study, the ERT method took approximately eight hours, the geomagnetic survey of the studied area took less than 20 minutes as part of a larger survey carried out in a 20 × 20 m grid that covered a broader extent of enclosure B and the GPR method required roughly the same time investment as the geomagnetic measurements.
In this work, we developed a methodology for analysing spatial correlations by selecting the main anomalies and using a binary system to assign a specific value to the points at which they are observed. As a result, we found that the ERT and GPR methods showed the best correlation, followed by the ERT and magnetic methods and, finally, the magnetic and GPR methods. These results were consistent with the values obtained for the Pearson correlation coefficient. We obtained numerical correlations at a particular depth, but data from different depths could be analysed in a similar way.
This study finds that spatial correlations are of fundamental importance for quantifying the relationships between the results obtained in multi-method studies and that correlation analyses must be carried out not only in a graphical way but also by applying a quantitative methodology.
Furthermore, we observe that the Pearson coefficient is an effective tool for investigating the correlations of geophysical anomalies, and it offers information about global correlation. To compare the spatial coincidence of the anomalies revealed by different methods, it is necessary to perform a quantitative correlation.
In addition, the results support the claim that the methodology that we have developed to quantify spatial correlations is suitable and clearly reflects the areas of coincidence of the anomalies identified using different methods.
Regarding the archaeological model of the site, the study area is part of a larger building with an N-S orientation in which we identified a central hall or distributor space of 15 × 6 m, surrounded by a series of rooms, some with internal subdivisions. The regular plan of the building, with a tripartite division centred on an axial corridor, offers close analogies with architectonic models of Roman/Italic origin. It fits very well with the early Roman chronology of the last occupation of the site, which is suspected to correspond to a military garrison during the Sertorian War (82-72 BC).
The ERT data indicate that the distributor space has a remarkable stratigraphic depth, detecting an anomaly of low resistivity that reveals a structure that reaches a depth of 2 m, which could be related to food storage and conservation. To the south, the area is delimited by the closing wall of the complex, probably built with slate basements and mud-brick walls, like the rest of the buildings already known. The magnetic and ERT data are highly consistent with the wall alignment present in other constructions that were found in excavated areas of the settlement. Regarding the roof, during previous archaeological surveys in this sector, several fragments of roof tiles and ceramic building material were recorded, including some tegulae raised flanges. These may be part of the areas of positive values obtained from the magnetic survey, corresponding to the roof collapse. Funding: This work has been partially supported by the research project "Desarrollo de métodos de mínima invasión para la revalorización socio-cultural de zonas arqueológicas" (PRI IB116150) funded by the Autonomous Government of Extremadura, by the project "Innovación abierta e inteligente en la EUROACE" with the reference 0049_INNOACE_4_E, by the European Union through the European Regional Development Fund included in COMPETE 2020, and through the Portuguese national funding agency for science, research and technology (FCT) in the scope of the projects SFRH/BSAB/143063/2018 and UIDB/04683/2020. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.