**1. Introduction**

Modern coastal landscapes are widely shaped by meandering fluvial, fluvio-tidal and tidal channels, which over the late Holocene accumulated complex and extensive sedimentary bodies. These bodies today often define subsoil permeable systems [1] that are often exploited as water reserves for agricultural, industrial and civil uses [2], and are extremely sensitive to saltwater intrusion [3,4] as well as to contamination [5,6]. These channelized bodies commonly consist of clean and poorly consolidated sand, which can also be a ffected by liquefaction processes [7]. The 2012 earthquake that occurred in the northeastern portion of the Po Plain (Italy) was the cause of sand eruptions that occurred along Holocene paleochannels and crevasse splay deposits down to an 8 m depth [8,9].

Aerial photographs and satellite images are excellent tools to identify and map late Holocene coastal channel networks since the former provide aerial data from the 1950s to present [10–13] and the latter provide multispectral analysis to highlight surficial paleochannel configurations [14–16].

Excellent examples of remote-sensing applications for these types of geomorphological studies can be found in the recent literature [17–21]. Despite the advantages, particularly in locating the position of these surficial bodies, and possibly distinguishing between tide- and fluvial-generated meanders [22], remote sensing alone is of course not capable of providing information at depth, thus remaining essentially a qualitative tool for the characterization of 3D geological structures. On the other hand, direct field surveys [23,24] and microrelief analysis [10,23,25] can provide further information either locally at depth or extensively at the surface. Regardless, a 3D reconstruction is still di fficult with these means only [26], if not as a result of interpolation of scarce scattered data.

Addressing these issues requires that extensive and high-resolution data are available to map large areas and, at the same time, investigate the subsoil to a certain depth. This calls for geophysical methods (as they are designed to collect data informative about the subsoil, unlike remote sensing sensu strictu) that can also be deployed rapidly with limited or no ground contact so that large areas can be investigated. The most suitable methods for these purposes are those based on electromagnetic processes. In particular, approaches based on electromagnetic induction (EMI) allow for noncontact subsurface investigation, with no intrinsic limitations as posed, for instance, to wave-propagation EM methods such as ground-penetrating radar (GPR) that can be strongly limited in their depth propagation by the ground electrical conductivity.

EMI is a well-established technique that dates back nearly one century [27] and is based on Faraday's law of electromagnetic induction. The technique is articulated in a variety of specific instrument designs and investigation strategies [28] ranging in investigation depth from very shallow (the first meter or so) to tens of kilometers. For shallow applications [29,30], EMI has had widespread use in hydrological and hydrogeological characterizations [31–33], hazardous waste characterization studies [34,35], precision-agriculture application locations [36–38], archaeological surveys [39,40], geotechnical investigations [41] and unexploded ordnance (UXO) detection [42]. EMI measurements at small scale are typically conducted in the frequency domain (frequency domain electromagnetics or FDEM), and the results are classically expressed as apparent electrical conductivities (ECa) [43] using the so-called low-induction number approximation [44]. In addition to ECa mapping, the development of multifrequency and multicoil instruments has recently enabled the possibility of inversion of EMI measurements to provide quantitative models of depth-dependent electrical conductivity (EC), as the di fferent acquisition configurations either in terms of coil geometry or frequency allow for multiple independent data to be acquired in su fficient number to warrant inversion. The majority of inversion algorithms use a 1D forward model based on either the linear cumulative sensitivity (CS) forward model proposed by [44] or nonlinear full solution (FS) forward models based on Maxwell's equations (e.g., [45,46]). As with EMI mapping, applications using inverted EMI data have also been diverse (e.g., [47–52]). Applications typically focus on using an inversion based on either the CS or an FS forward model to produce regularized, smoothly varying models of EC with fixed depths or sharply varying models of EC where layer depths are also a parameter. In the most advanced cases, a full 3D model of electrical conductivity can be reconstructed over a relatively large area, similar to what can be obtained at a larger scale by using, e.g., time-domain airborne EMI systems (e.g., [53]). The use of small FDEM measurement systems, with rapid response and easy integration into mobile platforms, is the key factor in the success of EMI techniques for near-surface investigations in these fields, as they allow dense surveying and real-time conductivity mapping over large areas in a cost-e ffective manner. However, su fficient control on the acquisition geometry is often needed, as the instrument response has a strong dependence also on the elevation above ground and the relative height of the primary and secondary coils [47].

The purpose of this paper is to show how the integrated use of remote sensing, EMI and direct stratigraphic investigations can provide an e ffective and comprehensive 3D view of the geometry of a

fluvial sedimentary body. Results from the present study highlight the importance of an integrated approach to understand subsurface deposits.

#### **2. Materials and Methods**

#### *2.1. Geological Setting and Study Area*

The Venetian Plain is located at the northeastern end of the Po Plain, the largest Italian alluvial plain, and was generated during Holocene transgression by aggradation of fluvial meandering channels [23,24]. Specifically, the study area is located at the boundary between the Venetian Plain and the Po River Delta, in a zone which is characterized by a dense network of alluvial ridges and sand bodies that are the geomorphological products of the complex interaction between the Adige and Po Rivers during the late Holocene (Figure 1a) [25,54]. These sedimentary bodies currently host a multilayered system of phreatic and confined aquifers that are a ffected by saltwater contamination [4,55] and intensive water exploitation. The fluvial sedimentation occurred in an aggrading setting related to the marine highstand, and meander belts often correspond to fluvial ridges slightly elevated over the floodplain (i.e., 2 to 5 m above sea level (asl)) [54]. The present surface is a typical lowland landscape, which developed in the last 5000 years by the avulsions of the Adige and Po rivers [25].

The investigated site is located near the village of Anguillara Veneta (Figure 1b), about 1 km north of the current channel of the Adige River, in an area with surface elevations ranging between 0.7 and 2.0 m asl, where traces of abandoned meanders are visible in several aerial images and could be followed for about 7 km, from Stanghella to Anguillara Veneta. These paleohydrographic traces run nearly parallel to the present Adige River, even if they are slightly out of the natural levee deposits connected to the fluvial ridge of Adige. The river activated its present direction since the early Middle Age, while before it used to flow along the meander belt running from Este to Monselice to Chioggia (from west to east) [56]. Near Anguillara Veneta, the present course of the Adige River cuts the so-called fluvial ridge of Rovigo–Saline–Cona, which was formed by the Po River between 4500 and 3500 years BP [25].

The area experienced strong anthropogenic activity since the Roman period, when extensive field systems were settled in the whole Venetian Plain and parts of the Po Delta [25]. A major phase of reclamation started in the 16th century, when the Venetian Republic started the strong managemen<sup>t</sup> of the river network, leading to the construction of the dense network of dikes, canals and ditches that still characterizes the landscape. During the same period, the Gorzone Canal was also cut, which represents the northern boundary of the study area, to convey the water discharge of the Agno–Guà–Frassine–Santa Caterina river system towards the Adriatic Sea. During the first part of the 20th century, the reclamation was extended to the coastal plain, where large portions of swamps and lagoon landscapes were drained through the excavation of canals and the use of pumping stations. These interventions made it possible to artificially lower the groundwater table below the surface and to cultivate seasonal crops (e.g., corn, wheat, bit and soya bean) and vegetables. In the last decades, several strong leveling interventions were carried out to improve the e fficiency of the new draining system, as in the field investigated for this study. Besides the positive results, unfortunately, the reclamation also induced fast land subsidence caused by groundwater withdrawal, compaction of the drained soil and degradation of the organic matter formerly present in the marsh sediments [57]. From the downlift rate of the natural subsidence, ranging around 1 mm/y [58], in the last century the velocity strongly increased up to average values of 2 to 5 mm/y with large sectors up to 10 mm/y [59].

A large number of historical photos and maps are available for the Venetian Plain, along with freely distributed satellite images. Study sites are easily accessible for geophysical investigation and sedimentary cores. For these reasons, this study area represents an ideal site to develop and test the proposed integrated approach.

**Figure 1.** The study site: (**a**) location of Anguillara Veneta town in the southern Venetian Plain, at the boundary with the Po Plain sensu strictu; (**b**) satellite image (2013) of the study area (yellow box) in Anguillara Veneta town, PD (Italy).

## *2.2. Remote Sensing*

The study site was selected after the identification of the paleomeanders in the aerial and satellite images. In particular, the first characterization was carried out on some satellite images available from Google Earth, which generally are pansharpened images of true-color composite bands of the Digital Globe company (Westminster, CO, USA) (e.g., Ikonos, QuickBird, WorldView and GeoEye missions). We selected the images with a detailed spatial resolution between 1.0–0.5 m (i.e., images 31/07/2004; 23/06/2007; 06/06/2010; 21/04/2012; 06/08/2013; 16/08/2013; 28/03/2015; 22/06/2017; 18/07/2018) and imported them into the GIS software ArcMap (version 10.7.1) [60] and QGIS 3.10 [61] for image processing and comparison with the images available as basemap reference in ESRI. Moreover, we also considered several zenithal conventional aerial pictures available from the cartographic service of Veneto Region [62], consisting of scanned versions of black/white and color pictures from 1955 to 2008, with scales from 1:33,000 to 1:5000.

To investigate the spectral characteristics of the field surface, we analyzed some images from the satellite Sentinel-2, obtaining the normalized di fference vegetation index (NDVI) and the normal di fference moisture index (NDMI) [63–65]. These indices have a geometric resolution of 10 m and, being sensitive to plant health and hydraulic stress, respectively [66,67], were used to improve the identification of the traces of paleomeanders by linking sedimentology to vegetation health of the area. In particular, we produced the NDVI and NDMI not only from summer scenes, but also from di fferent seasons and years by processing the multispectral bands through the semiautomatic classification plugin (SCP) [68] and raster calculator of QGIS 3.10. In fact, in the study area the cultivated plants change with seasons; in addition to the crops growing during summer, winter cultivations such as wheat, barley and di fferent vegetables can be present.

Satellite, aerial and processed images were visualised in a properly georeferenced 3D space provided by the Move 2018.2TM [69] software.

#### *2.3. Geophysical Investigations*

The EMI surveys at the site were collected using a GF Instruments CMD-Explorer probe [70] that is a six-coil system (three coplanar pairs) operating at a single frequency equal to 10 kHz. The probe can be operated in both horizontal coplanar (HMD) and vertical coplanar (VMD) configurations [71], providing six independent measurements that are generally associated with six di fferent apparent depths of investigation (Table 1). To acquire all six configurations for each geographical location it is, however, necessary to reoccupy with some acceptable degree of approximation the same location twice (once for each coil orientation).


**Table 1.** Technical specifications of the multicoil CMD Explorer FDEM.

The FDEM probe was mounted on a specifically designed wooden carriage and connected to a Trimble 5800 GPS for continuous positioning, collecting data every second [72]. The wooden support was towed by a small tractor (Figure 2a). The acquisition apparatus adopted satisfies two fundamental requirements that proved extremely e ffective in terms of data quality [73,74]:


**Figure 2.** Methodologies. (**<sup>a</sup>**–**b**) Geophysical acquisition: (**a**) the electromagnetic tool on the wooden sledge dragged by a small tractor on the study field and (**b**) the survey path. (**c**) Position of the recovered cores. Yellow dots and green triangles indicate the hand auger core and the drilled cores with rotary technique,respectively.

In this manner, we collected about 20,000 EMI data points (Figure 2b), each in both HMD and VMD modes, with about one point every 0.5 m along the acquisition lines. The lines have a mutual distance of roughly 10 m (Figure 2b).

EMI data were then inverted to retrieve real soil conductivity values. For this purpose, we used the Interpex IX1D inversion software [75], a 1D routine based on smooth depth inversion according to the so-called Occam's approach [76]. The very dense spatial sampling allowed for the reconstruction of the subsoil practically in a 3D fashion by the juxtaposition of the 1D inverted profiles. For all locations, the same number of layers (eight) was used for the inversion, thus producing a consistent dataset. The results are presented in terms of 2D horizontal maps at several depths, that were then georeferenced using the Move 2018.2TM [69] software, which also allowed creation of 3D surfaces, by the linear method, and tetravolumes.

## *2.4. Sedimentary Cores*

Six cores were recovered at the study site to analyze sedimentary features of the study deposits and provide ground-truth for geophysical and remote-sensing data (Figure 2c). The core locations were established based on remote-sensing and geophysical data. Cores were collected by using an Eijkelkamp hand auger and a continuous drilling core sampler with rotary technique. Three cores were recovered using an Eijkelkamp hand auger, through a gouge sampler with a length of 1 m and a diameter of 30 mm, which prevented sediment compaction. Depth for these cores spanned from 2.5 to 6 m (Figure 2c). Three additional cores were recovered using a continuous drilling core sampler with rotary technique. These latter cores, which were 10 cm in diameter and reached a depth of 10 m, were located in the upstream, central and downstream part of the study bar (Figure 2c). Cores were kept humid in PVC liners and successively cut longitudinally, sampled for grain-size analysis, photographed at high resolution and preserved for making dry peels with epoxy. Each core was characterized following the basic principles of facies analysis: highlighting sediment grain size, color, oxidation, sedimentary structures and occurrence of bioturbation, plant debris and shell fragments.

The terminology used in this work is graphically summarized in Figure 3. The channel thalweg is defined as the deepest part of the active channel, where the coarsest deposits occur. Ri ffles and pools are situated at bend inflections and bend apexes, respectively, and correspond to the shallower and the deeper portions of the thalweg, respectively. Sinuosity is calculated as the ratio between the along-channel distance between two adjacent ri ffles and their linear distance (Figure 3). Straight channels are characterized by sinuosity close to 1, whereas sinuous channels reach values higher than 2.5.

**Figure 3.** Descriptive terminology for fluvial meanders and related deposits: (**a**) alluvial plain morphological elements, (**b**) point-bar deposits and (**c**) meander-belt morphometry.
