Remote Sensing, Archaeological, and Geophysical Data to Study the Terramare Settlements: The Case Study of Fondo Paviani (Northern Italy)

During the Middle and Recent Bronze Age, the Po Plain and, more broadly Northern Italy were populated by the so-called “Terramare”, embanked settlements, surrounded by a moat. The buried remains of these archaeological settlements are characterized by the presence of a system of palaeo-environments and a consequent natural gradient in soil moisture content. These differences in the soil are often firstly detectable on the surface during the seasonal variations, with aerial, satellite, and Laser Imaging Detection and Ranging (LIDAR) images, without any information on the lateral and in-depth extension of the related buried structures. The variation in the moisture content of soils is directly related to their differences in electrical conductivity. Electrical resistivity tomography (ERT) and frequency domain electromagnetic (FDEM), also known as electromagnetic induction (EMI) measurements, provide non-direct measurements of electrical conductivity in the soils, helping in the reconstruction of the geometry of different buried structures. This study presents the results of the multidisciplinary approach adopted to the study of the Terramare settlement of Fondo Paviani in Northern Italy. Remote sensing and archaeological data, collected over about 10 years, combined with more recent ERT and FDEM measurements, contributed to the analysis of this particular, not yet wholly investigated, archaeological site. The results obtained by the integrated multidisciplinary study here adopted, provide new useful, interesting information for the archaeologists also suggesting future strategies for new studies still to be conducted around this important settlement.


Introduction
In the broad range of geophysical measurements applied to archaeology [1][2][3][4][5][6], electrical resistivity tomography (ERT) is today one of the most popular methods. Multichannel systems and new inversion techniques, developed between the 1980s and the 1990s, rapidly increased, in fact, the application and popularity of this technique in different archaeological contexts [7][8][9][10][11][12][13][14][15]. During the last 10 years, thanks to the development of 3D ERT surveys and new tools for 3D data inversion, this method enhances the possibility to reconstruct the spatial distribution and the shape of the archaeological targets, both in small or large areas, as well as in rural or urban context [16][17][18][19][20][21][22]. In this wide range of applications and possible uses of the ERT, undoubtedly, the main advantage offered by this method is  Plain has registered only a contraction of the population [51,52].
The buried remains of the Terramara, thanks to the presence of palaeo-environments and related natural differences in soil moisture content with the hosting system, represents an ideal target for ERT and FDEM measurements [53,54]. The first step in this recent multidisciplinary study of the Fondo Paviani settlement (Figure 2) moved from the pieces of evidence provided by the aerial photographs (Figures 2c and 3a-c).      The relative recent intensive agricultural practice, consisting of the regularization of the soil surfaces, in some cases, determines the removal of the upper part of the soil. This practice, also adopted at the Fondo Paviani in the last 40 years, could definitively compromise the archaeological deposit of the Terramare in some parts of the farm. The elevation data obtained by LIDAR ( Figure   3d) confirms local topographical modifications in the northern part of the site. Moving from these preliminary indications, a series of ERT lines collected in different parts of Fondo Paviani site help in  The comparison between an old aerial image ( Figure 3a) taken in 1955 and more recent images (Figure 3b,c) of the same area, respectively taken in 1990 and 2004, highlights the progressive attenuation of the northern boundaries of the settlement.
The relative recent intensive agricultural practice, consisting of the regularization of the soil surfaces, in some cases, determines the removal of the upper part of the soil. This practice, also adopted at the Fondo Paviani in the last 40 years, could definitively compromise the archaeological deposit of the Terramare in some parts of the farm. The elevation data obtained by LIDAR (Figure 3d) confirms local topographical modifications in the northern part of the site. Moving from these preliminary indications, a series of ERT lines collected in different parts of Fondo Paviani site help in the analysis of the actual condition and extension of these palaeo buried structures. The combination of ERT and FDEM measurements, also supported, in a test area, in the identification of the multitemporal, complex palaeo-systems spatial distribution. The integration between ERT and FDEM measurements, remote sensing data, and some pieces of evidence made by archaeological excavations, provided new interesting information about the studied area. The results of this multidisciplinary approach also support new hypotheses about the original structure of the Fondo Paviani embankment system, also suggesting the best way to continue the study around this important and not wholly investigated archaeological site.

Study Site Background
The broad fortified Terramara of Fondo Paviani, located in the Verona low plain, since 2007, is the focus of a multidisciplinary project [55] directed by the Prehistoric-Protohistoric research team at the Department of Cultural Heritage of the University of Padova. The settlement, developed between the last decades of the 14th and the 11th centuries BC, represents one of the most critical contexts for understanding the Late Bronze Age historical events of the Po Plain and, more broadly, in Northern Italy. Fondo Paviani was the central place of the complex political-territorial system known as Valli Grandi Veronesi polity. Between the 13th and the beginning of the 12th century BC, the settlement became the central keystone of the relations and trades connecting the Terramare Culture on the one hand with the Alpine area and with peninsular Italy, and on the other hand with Mycenean Greece, Cyprus, and the Levant [55].
Fondo Paviani also was one of the few centers that survived in the 12th century BC at the crisis producing the collapse of the Terramare Culture. Then it became the model for the development of the later population pattern and the new international trade network focused on the central place of Frattesina, on the Po river [56][57][58][59][60][61]. The development of the Fondo Paviani settlement reflects the typical model of the evolution of the Terramare-type settlements [62,63]. In the first phase, Fondo Paviani was surrounded only by a palisade and by a small ditch without remarkable water flows. During the 13th and the beginning of the 12th century BC a wide rampart was built and a ditch connected to the wetlands of the Menago river valley surrounded the site.
At that time, the settlement, which finally took the Terramara-type shape, reached the surface of 20 ha (16 ha of inner area and 4 ha of perimetral structures), an exceptional size for a Bronze Age settlement in Italy. The new fortified system, one of the biggest of the entire Terramare world, involved the active exploitation of pre-existing river channels and alluvial ridges, implying modifications and adaptations of the related landscape. Between the second half of the 12th and the 11th century BC, Fondo Paviani suffered a general crisis [60]. Broad areas of the settlement, for a long time occupied by houses and related infrastructures, were converted into cultivated spaces. In the last phase of the life of the settlement and many centuries after its abandonment, the wide rampart, although collapsed in several points and damaged by diggings, retained its primary function. The rampart represented an essential shield between the development of the stratification of the inner part of the settlement and the growing clogging of the ditch linking to the Menago valley [47].

Geomorphological Setting
The settlement of Fondo Paviani is part of the Valli Grandi Veronesi system, a lowland of Adige alluvial fan, on an alluvial ridge inside the valley of the Menago river belonging to the Bassa Pianura Veronese. The valley flows from the Adige alluvial fan through the Valli Grandi Veronesi, where it gets the shape of an N-W/S-E large depression with slight escarpments. The Menago valley is blocked by a large Late Pleistocene ridge ("Fabbrica dei Soci" ridge) immediately S-E of the site of Fondo Paviani ( Figure 2). This ridge, since the Holocene, has prevented the drainage of the valley, gradually changing the whole depression in a wetland, without covering the high bump of the settlement [64].
The Valli Grandi Veronesi area is well known [64][65][66][67][68] for the high preservation of the Bronze Age landscape. The depth of archaeological remains in this area is very shallow (between 0.2 and 1.4 m) [69,70]. Sometimes, the topography of the ancient structures corresponds to that of the current landscape [70], where this high preservation is the result of a lucky sequence of hydrogeological and anthropic events in the area from the end of Bronze Age to nowadays [65][66][67][68][69][70][71].
The main preservation factor of the Valli Grandi Veronesi ancient landscape is due to the presence of an alluvial/marshy clay layer [65], which gradually had covered the whole area from the Early Middle Age to the Modern Age. The presence of this layer and the general low drainage had made the Valli Grandi Veronesi a swamp almost uninhabited until the end of the 19th century when the reclamations of the area began.
For these reasons, the Valli Grandi Veronesi represents an optimal context both for archaeological research and for non-invasive analysis as aerial photo interpretation, DTM data from LIDAR, and geophysical measurements.
The hydrographic network contemporary to Fondo Paviani Terramara mainly consisted of spring-fed rivers flowing through the Menago Valley [65,66]. The sedimentary contribution of these rivers started between the end of the Early Bronze Age and the central phases of the Middle Bronze Age (between the 17th and the 16th century BC), is also at the origin of the formation of the ridge where, about two centuries later, the Fondo Paviani settlement developed [65]. A spring-fed river and Perteghelle channel, a watercourse certainly active between the Late Pleistocene and the beginning of the Holocene surround the N-W and S-W sides of the site. Probably, the Perteghelle palaeo-channel was not wholly inactive contemporarily to the life cycle of Fondo Paviani, although carrying slight flows of water [67].

Aerial Photograph Interpretation
The analysis and the consequent interpretation of aerial photographs represented the first step of the whole research. The aerial photograph interpretation helps in the identification of some hypothetical anthropic structures and natural features of the buried landscape, allowing in the plan of the field investigation the direct or non-direct verification of these features.  Table 1). The details of each aerial shot are reported in Table 1. Different years and months were considered to analyze the same area in different seasonal and cultivation conditions ( Table 1).
The high quality of the frames and the high visibility of the buried landscape reduce the processing to the enhancement and object classification steps [72,73]. The piece of evidence related to the archaeological and natural landscape features have been selected and compared, frame by frame, with different soil marks and crop marks [74]. A previous study of the Valli Grandi Veronesi landscape, based on aerial photo-interpretation [75], attributes light colors at features consisting of draining sediments (sand and coarse silt), often related to elevation structures, that can correspond to alluvial ridges or artificial ramparts, and dark colors at limited drainage areas and, more broadly, wetlands or fillings of abandoned water channels.
The aerial frames and the vector data of the photo interpretation were then loaded on Quantum GIS (release 2.18) and georeferenced. The contribution of GIS is essential to measure the features, to define their geometry better, and to determine their spatial relationships. Through GIS it was also possible to integrate the aerial photo interpretation with DTM data derived from LIDAR defining the geographical position of the features to plan further field investigations.

LIDAR
Thanks to the initiative of Consorzio di Bonifica Veronese, a LIDAR survey was carried out by C.G.R. (Compagnia Generale Riprese aeree) in Spring 2012. The survey was made from an elevation of 900 m using an airborne laser scanner equipped with Optech "Pegasus" sensors (elevation accuracy: 0.1 m; planimetric accuracy: 0.2 m; point density: 3 p/m 2 ). A DTM with a resolution of 0.5 m per cell was extracted from this dataset by C.G.R. These DTM results, therefore, already processed from LIDAR data, were used in this research. DTM data related to Fondo Paviani and its hinterland were then processed through GIS in order to better isolate the features recognized by aerial frame analysis. These data were also used to identify other elements of the ancient landscape that are not clearly evident through the photo interpretation (e.g., slight morphological discontinuities and slope breaks). At a later time, the analysis focused on the parts of the settlement fortified system that DTM can display, that is to say moat and rampart (remains of palisades cannot be readable through DTM). As for the moats, considering their negative shape, it is possible to have information from the DTM in terms of visibility and morphology, but not in terms of state of preservation. Instead, considering that ramparts are elevated structures that in larger Terramara settlements-as Fondo Paviani-could reach a width of 20 m and a height of 6 m, through DTM analysis it is possible to register not only data about presence and morphology, but also information about state of preservation. Moreover, DTM analysis has allowed to better plan the field investigations.

Stratigraphic Analysis
The stratigraphic analysis of the archaeological layers, carried out both in open area and in section, is one of the main topics of the Fondo Paviani project [55]. The open area investigations involved two excavation sectors in the North-Eastern part of the settlement: sector 1 (48 m 2 ), located between the inner part of the village and the rampart, and sector 2/2.1, positioned in the inner part of the village 38 m south-west from sector 1. Anthropogenic deposits, between 0.40 and 0.60 m deep from the ground level, were covered by the agrarian topsoil and by a sandy alluvial deposit dating to the first phases of the Iron Age. While sector 1 is still under investigation, the studies in sector 2/2.1, that ended in 2018, have allowed to reconstruct the living sequence, the housing structures, and other facilities, and more generally the activities, including handcrafting, carried out in this part of the settlement [55,61].
Current data from open area excavations do not allow to elaborate a complete planimetry of the dwellings and to add information on the fortified system, excepts in terms of chronological excursus; for this reason, open area investigations were not considered in this work.
The analysis of the E-W stratigraphic section (Figure 3), here discussed, and studied before in the frame of other projects [64], is fundamental in the comprehension of the fortified system of the Fondo Paviani settlement [66]. This archaeological section, 90 m long from a current drainage ditch crossing from E to W in a large part of the settlement (Figure 3), intercepts from W to E the inner stratification of the settlement, the perimetral fortification system, and a part of the wetlands of the Menago valley. The stratigraphic analysis of this section, carried out between 2007 and 2012 using a "genetic-process approach" [76], allowed to identify a complex sequence of natural and anthropogenic layers [55,66,77,78] covered by 0.4-0.6 m of agrarian topsoil. After this analysis, the stratigraphic sequence was measured and geographically positioned ( Figure 3). Pottery sherds and samples of organic material were collected from the section in order to date, both typologically and with 14C dating, the sequence, where many other samples supported micromorphological, archaeobotanical, and malacological analysis [55,66,78].

Electrical Resistivity Tomography (ERT)
In 2013 and 2015, an extensive field campaign of ERT acquisitions was planned at the archaeological site of Fondo Paviani to verify the perimeter of the settlement, its preservation, and its real extension. In total, five ERT sections were collected, the location of which, based on the pieces of evidence offered by remote data, is shown in the images in Figure 3. Three positions of interest were selected for the tomographies acquired in the SW-NE direction and two in the N-S direction ( Figure 3). The resistivity measurements were performed using an IRIS Syscal Pro 72 resistivity meter (Figure 4a). The acquisition dataset was optimized to take full advantage of the 10 physical channels available in this instrument. For the acquisition, a complete skip 4 (i.e., dipole spacing of five electrodes) dipole-dipole scheme was adopted, setting as a measurement criterion a pulse duration of 250 ms for each cycle and a target of 50 mV for potential readings for each current injection. In particular, L1, L6, L7 lines (  A quality factor "Q" (the difference between cycle results) equal to 5% was selected. For the ERT inversion, a regularized weighted least squares approach was adopted, according to the 'Occam's' rule [80] described in detail by LaBrecque et al. [81]. The smoothness of the resistivity distribution here calculated depends on the errors in the data set.
Binley et al. [82] demonstrated that a good evaluation of errors might be obtained by the analysis of these in direct and reciprocal measurements. Theoretically, these measurements shall be equal, providing the same resistance value. Any deviation may be interpreted as an error estimate, quantifying the error parameters useful for the inversion. In the present study, the inverted datasets present errors smaller than 5%. The visualization of the inverted data was done using Surfer 10 software (Golden Software). A quality factor "Q" (the difference between cycle results) equal to 5% was selected. For the ERT inversion, a regularized weighted least squares approach was adopted, according to the 'Occam's' rule [80] described in detail by LaBrecque et al. [81]. The smoothness of the resistivity distribution here calculated depends on the errors in the data set.
Binley et al. [82] demonstrated that a good evaluation of errors might be obtained by the analysis of these in direct and reciprocal measurements. Theoretically, these measurements shall be equal, providing the same resistance value. Any deviation may be interpreted as an error estimate, quantifying the error parameters useful for the inversion. In the present study, the inverted datasets present errors smaller than 5%. The visualization of the inverted data was done using Surfer 10 software (Golden Software).

Frequency Domain Electromagnetic (FDEM)
In the frame of the geophysical measurements carried out over a few years at Fondo Paviani, in 2015 a small area was selected with clear traces of potential structures belonging to the perimeter system of the embanked site, visible from aerial and satellite images (Figure 3).
In this area, an extensive FDEM multifrequency measurement was made and integrated by the L6 ERT line acquisition, as shown in Figure 3.
For FDEM data collection, a Geophex GEM-2 multifrequency conductivity meter was used ( Figure 4b). The instrument operates with a fixed source-receiver separation of 1.66 m at multiple frequencies ranging from 330 Hz to 48 kHz. Both quadrature and in-phase responses were recorded carrying the instrument in the vertical-dipole configuration at a height of 1 m above the ground surface and using seven frequencies: f 1 = 775 Hz, f 2 = 3775 Hz, f 3 = 6275 Hz, f 4 = 10,475 Hz, f 5 = 17,275 Hz, f 6 = 29,025 Hz, and f 7 = 47,025 Hz. The FDEM data were collected with the GEM-2 ski oriented in-line to the walking path, at about 0.5-m sampling interval along with parallel profiles, approximately 2 m apart, in a northwest-southeast direction. Synchronizing the GEM-2 device with a differential GPS receiver (Trimble 5800), UTM coordinates of each FDEM measurement point were recorded with sub-meter accuracy. Raw data were preliminarily analyzed for detecting and removing possible DC (static) shifts, outliers, and short-wavelength noises, which usually adversely affect the quality of the inversions, causing spiking and very different solutions with adjacent soundings. Since the in-phase component of the data revealed too noisy, the only quadrature component of the filtered data was inverted to produce a 1D electrical conductivity profile below each measurement point. Then, the 1D models obtained were stitched together to build a pseudo-3D volume of the investigated area. Inversion was performed using the FDEMtools [83], a free MATLAB software package implementing the numerical algorithms mainly discussed by Deidda et al. [84][85][86]. A layered starting model consisting of 30 layers, to a depth of 3.5 m, was used to invert the electromagnetic data. During the inversion, the magnetic permeability of the layers, as well as the depth of their interfaces, remain fixed.

The Nonlinear Forward Problem and the Inversion Procedure
The forward modeling used to calculate the nonlinear EM response of a layered half-space ( Figure 5) for dipole source excitation is well known [87][88][89]. It is based on Maxwell's equations, suitably simplified thanks to the cylindrical symmetry of the problem, since the magnetic field sensed by the receiver coil is independent of the rotation of the instrument around the vertical axis. When the axes of the coils are aligned vertically to the ground, the model prediction M(σ, µ), defined as the ratio of the secondary (H S ) to the primary (H P ) EM field, can be expressed in terms of Hankel transforms of order zero: where λ is a variable of integration with no particular physical meaning, h is the height of the instrument above the ground, r the coil separation, J 0 the Bessel function of order 0, and R 0 (λ) the response kernel, which is a complex value function of the parameters that describe the layered subsurface (i.e., for the k-th layer: the electrical conductivity σ k , the magnetic permeability µ k , and the layer thickness d k ) besides the frequency and λ.
where λ is a variable of integration with no particular physical meaning, h is the height of the instrument above the ground, r the coil separation, J0 the Bessel function of order 0, and R0(λ) the response kernel, which is a complex value function of the parameters that describe the layered subsurface (i.e., for the k-th layer: the electrical conductivity σk, the magnetic permeability μk, and the layer thickness dk) besides the frequency and λ. The nonlinear inversion procedure proposed by Deidda et al. [84][85][86] is a general procedure that allows the estimation of the electrical properties (electrical conductivity and magnetic permeability) of the subsurface by inverting the complex multi-depth response of different electromagnetic devices designed to record data at multiple coil spacings, using a single frequency, or at multiple frequencies with a fixed coil spacing. Let us suppose that the aim of the survey is estimating the vertical distribution of the electrical conductivity, using a multifrequency electromagnetic dataset (e.g., data recorded with the GEM-2 device). Fixing the magnetic permeability µ k (k = 1, . . . , M), the best approximation of conductivity values σ k (k = 1, . . . , M) can be found minimizing the Euclidean norm of the complex residual r(σ) between the data b(ω i ) (where ω i with i = 1, . . . , N, are the operating angular frequencies) and the forward model prediction based on Equation (1), that is, Alternatively, either the in-phase (real) or the quadrature (imaginary) component of the residual r(σ) could also be minimized. The nonlinear minimization problem is solved with an algorithm based on a damped regularized Gauss-Newton method. Assuming the Fréchet differentiability of r(σ), the problem is linearized at each iteration by means of a first order Taylor expansion with the analytical Jacobian (sensitivity matrix), which makes the computation faster and more accurate than using a finite difference approximation [84,85]. The regularized solution to each linear subproblem is then computed by the truncated generalized singular value decomposition (TGSVD) [90], employing different regularization operators (first and second derivatives).

Results
Aerial photo interpretation allowed us to preliminary define both anthropic and natural elements that probably delimited the settlement area of Fondo Paviani. Stratigraphic analysis of the archaeological section confirmed that the light and straight mark, visible from aerial photos in the N-E part of the settlement, correspond to the structure of an impressive artificial rampart. The dark areas E of the rampart correspond to the Menago valley wetlands ( Figure 6). Even the elevation model (DEM) obtained from LIDAR confirms the coherence of these data with the presence of both the rampart and the Menago valley depression. However, only the stratigraphic analysis of the archaeological section allowed us to fully understand the characteristics of the fortified system and its evolution during the centuries in the first phase, which started in the last decades of the 14th century BC [47], the settlement was not enclosed by a rampart, but probably only by a palisade and a small ditch (Figure 6b) [55]. The wide rampart was built between the 13th and the beginning of the 12th century BC [47]. The structure of this one consisted of a silt-sandy basal core and a top body of fine-textured sediments, built using an imposing internal wooden frame (Figure 6b) [55]. Thanks to the aerial photo interpretation, an alternating sequence of light and dark polylinear marks have been identified along the S-E and S-W sides of the settlement. The outline and shape of these marks are compatible both with anthropic structures, such as ramparts and ditches, and with natural morphologies, such as alluvial ridges and channels. The absence of recent ground Even the elevation model (DEM) obtained from LIDAR confirms the coherence of these data with the presence of both the rampart and the Menago valley depression. However, only the stratigraphic analysis of the archaeological section allowed us to fully understand the characteristics of the fortified system and its evolution during the centuries in the first phase, which started in the last decades of the 14th century BC [47], the settlement was not enclosed by a rampart, but probably only by a palisade and a small ditch (Figure 6b) [55]. The wide rampart was built between the 13th and the beginning of the 12th century BC [47]. The structure of this one consisted of a silt-sandy basal core and a top body of fine-textured sediments, built using an imposing internal wooden frame (Figure 6b) [55]. Thanks to the aerial photo interpretation, an alternating sequence of light and dark polylinear marks have been identified along the S-E and S-W sides of the settlement. The outline and shape of these marks are compatible both with anthropic structures, such as ramparts and ditches, and with natural morphologies, such as alluvial ridges and channels. The absence of recent ground checks in these sectors of the settlement prevented from establishing the effective origin (natural or artificial) of these marks. However, the internal light marks, having dimensions comparable to the traces of the N-E rampart, can correspond to the remains of an artificial embankment. The N-W side of the settlement seems to be protected by the natural alluvial ridge embankment offered by the so-called Perteghelle palaeo-channel, visible in the aerial frames, and highlighted by the DTM extracted from LIDAR data.
The twisting linear mark, visible both by aerial photo interpretation and DTM analysis and that cross the settlement from North to South, can be interpreted as a palaeo-channel. Focused stratigraphic analysis has allowed to date this channel to Middle Ages.
The results of the ERT measurements are shown in the images in Figure 7.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 checks in these sectors of the settlement prevented from establishing the effective origin (natural or artificial) of these marks. However, the internal light marks, having dimensions comparable to the traces of the N-E rampart, can correspond to the remains of an artificial embankment. The N-W side of the settlement seems to be protected by the natural alluvial ridge embankment offered by the socalled Perteghelle palaeo-channel, visible in the aerial frames, and highlighted by the DTM extracted from LIDAR data. The twisting linear mark, visible both by aerial photo interpretation and DTM analysis and that cross the settlement from North to South, can be interpreted as a palaeo-channel. Focused stratigraphic analysis has allowed to date this channel to Middle Ages.
The results of the ERT measurements are shown in the images in Figure 7. To allow smooth and fast comparison between the different tomographies, the resistivity range of the ERT sections, acquired in different seasons and years, was standardized between 0 and 50 Ohm*m. The ERT sections in Figure 7 show the resistivity distribution in five areas along the perimeter of the Fondo Paviani settlement. The ERT results highlight the current state of the studied embanked system and its high heterogeneity, both lateral and in-depth, from area to area. L1 was collected corresponding to the position of the documented stratigraphic section. In Figure 8, the higher resistivity anomaly (red), localized in the NE part of the L1 section, identifies the preserved rampart. The L1 section was then used as a marker or term of comparison to evaluate, using other ERT sections, the degree of preservation of the boundary system in different areas. In general, ERT To allow smooth and fast comparison between the different tomographies, the resistivity range of the ERT sections, acquired in different seasons and years, was standardized between 0 and 50 Ohm*m. The ERT sections in Figure 7 show the resistivity distribution in five areas along the perimeter of the Fondo Paviani settlement. The ERT results highlight the current state of the studied embanked system and its high heterogeneity, both lateral and in-depth, from area to area. L1 was collected corresponding to the position of the documented stratigraphic section. In Figure 8, the higher resistivity anomaly (red), localized in the NE part of the L1 section, identifies the preserved rampart. The L1 section was then used as a marker or term of comparison to evaluate, using other ERT sections, the degree of preservation of the boundary system in different areas. In general, ERT L1 (Figures 7 and 8) identifies a transition from more resistive to low resistive values moving from SW to NE in the last part of the section. In particular, a clear conductive deep incision is here registered in the NE part, below the identified high resistive anomaly corresponding to the preserved rampart. The ERT L2-3 section (Figure 7), also collected in the E sector, south of the L1 line and almost parallel to it, shows a low lateral variation of the resistivity compared to L1 with a shallow conductive part. The two other ERT sections L6 and L7 were acquired parallel to each other (Figure 7). These describe a general coherent system, with a more resistive area in the northern part, and the distribution of shallow conductive bodies in the southern part of the section.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 22 The archaeological site of Fondo Paviani, characterized by the presence of the remains of a housing settlement consisting of a perimeter system, with a moat, characterized by river deposits, and a bright contrast to the hosting system, represents in this sense the ideal case-study for these types of non-invasive measurements. Thanks to the information obtained from the L1 section ( Figures  7 and 8), with the identification of the response of the preserved rampart, detected by the stratigraphic section made in the same position (Figure 8), it was, in fact, possible to verify what was hypothesized by the analyses carried out on aerial images and by the LIDAR data. The ERT sections, mainly made in the northern sector of the embanked settlement, confirm the absence of the preserved rampart in this sector in the shallow subsoil, where one would expect to find the same, as in section L1. Therefore, it can be reasonably said that in the sections corresponding to the ERT lines from L2-3 to L7 the highest and therefore most shallow system relating to the rampart has been removed due to the agricultural practices that have taken place in the last 40 years. This data would also be confirmed by the topographical data extracted by LIDAR (Figure 3), where mean lower altitudes are registered corresponding to the areas in which the ERT sections were acquired. However, from the comparison between the DTM data extracted from LIDAR and the ERT sections, it is possible to make further considerations on the possible correlation between the topographical surface data and the real buried deep geomorphological data. In particular, despite the apparent identical topographic information relating to the NE sector at the two ERT sections L1 and L2-3, parallel to each other, the geophysical data shows that the L1 section intercepts a deep conductive incision where instead the section L2-3 intersects a shallower system. From these real data, it is, therefore, possible to affirm that the topographical data and the overall view, albeit informative of the DTM data, do not always correspond to the real deep-buried structure. The ERT sections undoubtedly provided additional information of extreme interest regarding the geomorphology of the area. In particular, only the lines L1 and L7 seem to intercept some deep incisions characterized by conductive deposits for the whole investigation depth. In contrast, in the other ERT sections, the conductive structures, attributable to the presence of river deposits related to the perimeter system of the dam site, appear less in-depth.
In the areas of the settlement not investigated by geophysical measurements, the DTM analysis allowed us to monitor the state of preservation of the elevated structures (natural or anthropogenic) Finally, ERT section L4-5 (Figure 7), collected in the supposed western boundary of the settlement, describes a low lateral variation of the resistivity, similar to this registered in the section L2-3, except to an apparent resistive deep incision. No evidence about the presence of a preserved rampart are registered in all ERT lines from L2-3 to L6. The inversion of FDEM data made it possible to extract a depth map of the distribution of electrical conductivity 2 m below the surface in the area where the ERT L6 was acquired. Figure 9 shows the location of the investigated area with the FDEM method and the relative position of the ERT L6 line. The electrical conductivity data of the FDEM map was then converted into resistivity values to allow us a direct comparison with the ERT L6 line (Figure 9), using its real data range (0-100 Ohm*m) to evaluate the degree of information of these two different methods and their consistency.
The resistivity pattern visible in the FDEM map obtained by the inversion of this dataset, highlights and confirms the presence of the palaeo-system visible from the satellite image (Figure 9a). The ERT L6 section, entirely consistent with the FDEM data, completes with more details the information provided by the EM method about lateral extension and total depth of these buried geomorphological palaeo-structures. measurements carried out in a test area, in which the L6 ERT section was acquired, thanks to the inversion of the data, demonstrate absolute consistency with the ERT measurements and, therefore, the possibility of integrating this different kind of data. ERT combined with inverted FDEM data, thus confirming the possibility to quickly and precisely verify the extension of the structures only locally identified with ERT ( Figure 9). The integration between ERT and FDEM methods allows us to also answer question c positively. Undoubtedly, the two geophysical methods allow us to detect the different targets in the subsoil and, therefore, to define the real stratigraphic differences not detectable from the aerial photo images, where the surface signal does not precisely define the nature and the stratigraphic relationship between the different targets.
Finally, the integrated multidisciplinary approach adopted in this study, answers to point (d), highlighting how the information provided by every single method used, although fundamental, may be limited in the absence of a multidisciplinary approach as adopted here. The positive evidence obtained thanks to this type of integration allows us to outline future investigation strategies to be adopted to complete the study of this important embanked site.

Conclusions
The multidisciplinary study carried out on the Terramara of Fondo Paviani, with the contribution of remote sensing data, geophysical data, and the analysis of a stratigraphic section, provides new interesting information about this significant settlement of the Late Bronze Age. The first advantage of the combination of remote sensing data with geophysical data is the possibility of reconstructing the real geometry of the buried structures and their spatial distribution and the consequent stratigraphic relationship between them. It is clear that the evidence of the presence of the palaeo-channels, readable with remote sensing data, loses the information about the stratigraphy

Discussion
The archaeological questions underlying the multidisciplinary study conducted in the Bronze Age embanked site of Fondo Paviani and described in this contribution can be summarized in the following points: (a) Check the state of preservation of the perimetric system in different parts of the site; (b) Verify the degree of information available by the combination of different geophysical methods measuring the differences in the electrical conductivity of the soils and therefore defining the geometry (size and depth) and the spatial distribution of the palaeo-channels that characterize this embanked system; (c) Identify the stratigraphic/sedimentation relationships describing the reactivation sequence of older river systems with more recent ones; (d) Define a multidisciplinary protocol applicable to the whole site, capable of providing clear information on the position of the structures of interest.
To answer the question of point (a), we started from the evidence made by the analysis of some aerial photographs relating to a chronological span of about 40 years, which have also guided the excavations and archaeological reconnaissance for over a decade. Over time, aerial photographs have shown a marked decrease in the contrast, which was previously visible between the traces of paleochannels in the northern perimeter of the site and the hosting system. The stratigraphic section drawn up following the excavation carried out in a part of the perimeter of the system, which is preserved today, constituted the first direct useful data to plan the field strategy for the geophysical measurements carried out with the use of the ERT and the FDEM methods.
The ERT sections acquired in different points of the embankment system of Fondo Paviani allow here to formulate some interesting considerations, made possible thanks to the comparison and calibration of the no-direct data with the direct ones provided by the stratigraphic section. In general, ERT measurements find a wide application and allows us to obtain excellent results in contexts and on systems where there is a marked gradient of the soil water content between the target and the hosting system.
The archaeological site of Fondo Paviani, characterized by the presence of the remains of a housing settlement consisting of a perimeter system, with a moat, characterized by river deposits, and a bright contrast to the hosting system, represents in this sense the ideal case-study for these types of non-invasive measurements. Thanks to the information obtained from the L1 section (Figures 7  and 8), with the identification of the response of the preserved rampart, detected by the stratigraphic section made in the same position (Figure 8), it was, in fact, possible to verify what was hypothesized by the analyses carried out on aerial images and by the LIDAR data. The ERT sections, mainly made in the northern sector of the embanked settlement, confirm the absence of the preserved rampart in this sector in the shallow subsoil, where one would expect to find the same, as in section L1.
Therefore, it can be reasonably said that in the sections corresponding to the ERT lines from L2-3 to L7 the highest and therefore most shallow system relating to the rampart has been removed due to the agricultural practices that have taken place in the last 40 years. This data would also be confirmed by the topographical data extracted by LIDAR (Figure 3), where mean lower altitudes are registered corresponding to the areas in which the ERT sections were acquired. However, from the comparison between the DTM data extracted from LIDAR and the ERT sections, it is possible to make further considerations on the possible correlation between the topographical surface data and the real buried deep geomorphological data. In particular, despite the apparent identical topographic information relating to the NE sector at the two ERT sections L1 and L2-3, parallel to each other, the geophysical data shows that the L1 section intercepts a deep conductive incision where instead the section L2-3 intersects a shallower system. From these real data, it is, therefore, possible to affirm that the topographical data and the overall view, albeit informative of the DTM data, do not always correspond to the real deep-buried structure. The ERT sections undoubtedly provided additional information of extreme interest regarding the geomorphology of the area. In particular, only the lines L1 and L7 seem to intercept some deep incisions characterized by conductive deposits for the whole investigation depth. In contrast, in the other ERT sections, the conductive structures, attributable to the presence of river deposits related to the perimeter system of the dam site, appear less in-depth.
In the areas of the settlement not investigated by geophysical measurements, the DTM analysis allowed us to monitor the state of preservation of the elevated structures (natural or anthropogenic) on a preliminary basis. In particular, in north-eastern and south-western corners of the settlement, an abrupt and regular lowering of the surfaces is visible. These anomalies, considering that they follow exactly agricultural field borders, can be interpreted as the outcome of a strong agricultural impact with possible partial destruction of the ancient landscape features.
Considering the point b), the ERT method, capable of providing data on the lateral and in-depth extension of the targets, could be integrated into vast areas by faster soil electrical conductivity measurements. These can be obtained thanks to the use of multi-frequency FDEM systems that allow the inversion of the data providing detailed information at different depths. Multi-frequency FDEM measurements carried out in a test area, in which the L6 ERT section was acquired, thanks to the inversion of the data, demonstrate absolute consistency with the ERT measurements and, therefore, the possibility of integrating this different kind of data. ERT combined with inverted FDEM data, thus confirming the possibility to quickly and precisely verify the extension of the structures only locally identified with ERT ( Figure 9).
The integration between ERT and FDEM methods allows us to also answer question c positively. Undoubtedly, the two geophysical methods allow us to detect the different targets in the subsoil and, therefore, to define the real stratigraphic differences not detectable from the aerial photo images, where the surface signal does not precisely define the nature and the stratigraphic relationship between the different targets.
Finally, the integrated multidisciplinary approach adopted in this study, answers to point (d), highlighting how the information provided by every single method used, although fundamental, may be limited in the absence of a multidisciplinary approach as adopted here. The positive evidence obtained thanks to this type of integration allows us to outline future investigation strategies to be adopted to complete the study of this important embanked site.

Conclusions
The multidisciplinary study carried out on the Terramara of Fondo Paviani, with the contribution of remote sensing data, geophysical data, and the analysis of a stratigraphic section, provides new interesting information about this significant settlement of the Late Bronze Age. The first advantage of the combination of remote sensing data with geophysical data is the possibility of reconstructing the real geometry of the buried structures and their spatial distribution and the consequent stratigraphic relationship between them. It is clear that the evidence of the presence of the palaeo-channels, readable with remote sensing data, loses the information about the stratigraphy and, therefore, historical information about the time sequence of these structures. The singular information provided from remote sensing data, especially in the archaeological context, such as the one analyzed here, inserted in an active agricultural system, also does not allow us to define the real state of preservation of the buried structures. Additionally, the remote sensing data clarifies the development dynamics of the settlement even less, based on the exploitation of natural waterways and partly on the construction of artificial structures, that can be clarified only with the contribution of archaeological stratigraphic analysis and partially with geophysical data. The results of the multidisciplinary study carried out at the Terramara of Fondo Paviani, here presented, therefore, allows us to obtain important information on the state of preservation of this settlement, also suggesting some preliminary considerations on its fortification structure and evolution, and proposing new perspectives and new strategies for future research on this important archaeological site.