*2.3. 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/m2). 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.

## *2.4. 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 m2), 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].

#### *2.5. 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 o ffered 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 (Figure 3) were acquired using 48 electrodes with 2 m spacing, for a total length of 94 m and an investigation depth of approximately 18 m. Lines L2-3 and L4-5 (Figure 3) are instead the results of two acquisitions with 72 electrodes, spacing 2 m, for a total length of 142 m, and an investigation depth of about 28 m. For each ERT line, direct and reciprocal measurements were registered by swapping current with potential electrodes, to estimate the errors in the dataset [79].

**Figure 4.** (**a**) ERT instrumentation during an acquisition at Fondo Paviani site; (**b**) Multifrequency conductivitymeter used for the measurements at Fondo Paviani site.

A quality factor "Q" (the di fference 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).

#### *2.6. 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: f1 = 775 Hz, f2 = 3775 Hz, f3 = 6275 Hz, f4 = 10,475 Hz, f5 = 17,275 Hz, f6 = 29,025 Hz, and f7 = 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–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–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(<sup>σ</sup>, <sup>μ</sup>), defined as the ratio of the secondary (*HS*) to the primary (*HP*) EM field, can be expressed in terms of Hankel transforms of order zero: 

$$\frac{H\_S}{H\_P} = \int\_0^\infty \lambda^2 e^{-2h\lambda} \mathcal{R}\_0(\lambda) J\_0(r\lambda) d\lambda,\tag{1}$$

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 λ.

**Figure 5.** Schematic representation of the subsoil discretization and parameterization along with the coils of the measuring device above the ground.

The nonlinear inversion procedure proposed by Deidda et al. [84–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**(<sup>ω</sup>*i*) (where ω*i* with *i* = 1, ... , *N*, are the operating angular frequencies) and the forward model prediction based on Equation (1), that is,

$$\widehat{\boldsymbol{\sigma}} = \arg\min\_{\boldsymbol{\sigma} \in \mathbb{R}^M} \frac{1}{2} \left\| \mathbf{r}(\boldsymbol{\sigma}) \right\|^2 = \arg\min\_{\boldsymbol{\sigma} \in \mathbb{R}^M} \frac{1}{2} \left\| \mathbf{b}(\boldsymbol{\omega}\_i) - \boldsymbol{\mathcal{M}}(\boldsymbol{\sigma}) \right\|^2. \tag{2}$$

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 **<sup>r</sup>**(σ), 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).
