Geophysical and Sedimentological Investigations Integrate Remote-Sensing Data to Depict Geometry of Fluvial Sedimentary Bodies: An Example from Holocene Point-Bar Deposits of the Venetian Plain (Italy)

Abstract

Over the past few millennia, meandering fluvial channels drained coastal landscapes accumulating sedimentary successions that today are permeable pathways. Propagation of pollutants, agricultural exploitation and sand liquefaction are the main processes of environmental interest affecting these sedimentary bodies. The characterization of these bodies is thus of utmost general interest. In this study, we particularly highlight the contribution of noninvasive (remote and ground-based) investigation techniques, and the case study focuses on a late Holocene meander bend of the southern Venetian Plain (Northeast Italy). Electromagnetic induction (EMI) investigations, conducted with great care in terms of sonde stability and positioning, allowed the reconstruction of the electrical conductivity 3D structure of the shallow subsurface, revealing that the paleochannel ranges in depth between 0.8 and 5.4 m, and defines an almost 260 m-wide point bar. The electrical conductivity maps derived from EMI at different depths define an arcuate morphology indicating that bar accretion started from an already sinuous channel. Sedimentary cores ensure local ground-truth and help define the evolution of the channel bend. This paper shows that the combination of well-conceived and carefully performed inverted geophysical surveys, remote sensing and direct investigations provides evidence of the evolution of recent shallow sedimentary structures with unprecedented detail.

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 affected 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,11,12,13] and the latter provide multispectral analysis to highlight surficial paleochannel configurations [14,15,16].

Excellent examples of remote-sensing applications for these types of geomorphological studies can be found in the recent literature [17,18,19,20,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 difficult 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,32,33], hazardous waste characterization studies [34,35], precision-agriculture application locations [36,37,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 different acquisition configurations either in terms of coil geometry or frequency allow for multiple independent data to be acquired in sufficient 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,48,49,50,51,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-effective manner. However, sufficient 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 effective 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 affected 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 management 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 efficiency 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.

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 difference vegetation index (NDVI) and the normal difference moisture index (NDMI) [63,64,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 different 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 different vegetables can be present.

Satellite, aerial and processed images were visualised in a properly georeferenced 3D space provided by the Move 2018.2^TM [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 different apparent depths of investigation (

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

Instrument Probe	Coil Interdistance (m)	Frequency	Nominal Exploration Depth (Horizontal Mode HMD /Vertical Mode VMD)
1	1.48	10 kHz	2.2/1.1 m
2	2.82	10 kHz	4.2/2.1 m
3	4.49	10 kHz	6.7/3.3 m

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

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 effective in terms of data quality [73,74]:

(a)

Reoccupation of the same location is warranted by the GPS within the required precision (note that the sonde is a few meters long);

(b)

The setting of the sonde is the same at all locations, with no changes of either the height above ground or the setting of the sonde that is maintained largely horizontal.

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.2^TM [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. Riffles 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 riffles and their linear distance (Figure 3). Straight channels are characterized by sinuosity close to 1, whereas sinuous channels reach values higher than 2.5.

3. Results

3.1. Remote Sensing

The analyzed satellite images and aerial pictures give a consistent idea of the changes in the land use of the area and the visibility of the paleohydrographic traces over the last 60 years. In particular, recent high-resolution true-color composite images allow mapping these features with submetric accuracy. The traces mainly consist of cropmarks with vegetation suffering from hydraulic stress during the growing season, and dead patches in July and August. Significant changes can be observed over very short periods (compare, e.g., Figure 4d,f) during the hot season, while sensitivity is low in spring time (e.g., Figure 4e) and is lost in winter, when the fields are covered by scarce vegetation and are bare and plowed, respectively. By comparing satellite images mostly taken during the growing season (i.e., those that best show the traces) differences in vegetation colors allow the identification of buried morphologies of two distinct reaches, hereafter named Reach A and Reach B, and a crevasse splay. Reach A consists of two major bends, called Bend 1 and Bend 2 (Figure 4b). The NDVI and NDMI images derived from Sentinel-2 summer images (Figure 5c,d) clearly show the differential behavior of the paleomeanders in comparison to the external floodplain and the inner portion surrounded by Reach A; the vegetation growing on more permeable sandy soil is less healthy than the one living on finer sediments. The paleohydrographic traces are much more evidenced by NDVI and NDMI during the summer season and they clearly help in identifying the general pattern of the abandoned channels. However, the rather low resolution of the Sentinel-2 images does not contribute significantly to discriminating in detail the specific morphological and sedimentological features composing the paleohydrographic traces.

Bend 1 is an open bend, with an SSW–NNE bend axis, and is characterized by a sinuosity of about 2.2 and a radius of curvature of ca. 140 m. The scroll-bar pattern is particularly clear from cropmarks, in the northcentral portion of the bend (Figure 4c,d), showing a different signal compared to the residual channel fill (i.e., light cropmark when the others are dark, and vice-versa), testifying the progressive growth of the meander bend. The channel fill displays a width of about 15–20 m and can be better defined where bounded by opposite-trending scroll-bar patterns, like in the upstream side of the bend. The riffle-to-riffle distance on the channel fill is about 260 m.

Bend 2 is a low sinuosity bend (i.e., 1.12), with an NW–SE bend axis, characterized by an estimated radius of curvature and riffle-to-riffle distance of 135 and 230 m, respectively. Bend 2 is sited upstream from Bend 1, and shows a scroll-bar pattern that testifies a progressive expansional growth style (sensu [77]) of the bend.

Reach B forms a bend occurring south of Reach A, but is less visible from satellite images and its sinuosity cannot be defined. The radius of curvature is ca. 350 m and the axis of the bend trends ca. SW–NE, although satellite images do not show a clear bar-scroll pattern and the position of the relevant channel fill. Reach B cuts over Reach A suggesting that it developed after a chute channel that cut off Reach A [78], which was later abandoned. Additionally, along the eastern side of Reach B, (Figure 4d,f) a divergent pattern of minor channels point to a local development of a crevasse splay sourced from the downstream side of the bend. Several straight dark stripes, with a width of about 1 m, located on Reach B, can be interpreted as traces of abandoned ditches that were associated with a drainage system dating back to Renaissance times and dismissed later on (Figure 4d).

3.2. 3D Electrical Conductivity Model from EMI

The inversion of the EMI data produced a 3D volume of electrical conductivity values. The results are shown in Figure 6, where we elected to show the volume sliced horizontally at eight depth levels down to a depth of nearly 8 m below ground, corresponding each to a layer selected in the inversion approach. Note that the inversion was conducted with an Occam approach, but using a limited number of layers compatible with the information contained in the six different acquisition configurations obtainable with the CMD Explorer instrument.

An arcuate sedimentary body having low electrical conductivity (i.e., a resistive body) is clearly visible at a depth between 1 and 6 m below ground. The internal boundary of this arcuate body (see Slice 5 in Figure 6b) is fully visible in the maps, and shows a radius of curvature and a sinuosity of ca. 60 and 2.3, respectively. The external boundary (see Slice 5 in Figure 6b) slightly debouches from the maps, but its radius of curvature and sinuosity can be estimated to be ca. 135 and 2.2, respectively (Figure 6b), as also confirmed by remote-sensing results. Orientation of the outer boundary of this body fits with the orientation of meander Bend 1 of Reach A, as depicted by remote-sensing analyses, and is also consistent with the associated scroll pattern (Figure 4b), suggesting that these low-resistivity deposits represent the point-bar body associated with meander Bend 1 of Reach A. Of course, the main contribution of the EMI data is to provide continuous and extensive depth information that is not available from remote sensing. In the shallower layers (Slices 1–5), the arcuate point-bar body presents low conductivity values with σ < 20 mS/m, and its conductivity is still close to 40 mS/m at about 5–6 m below ground (Slices 6 and 7; Figure 6b). Note that the width of the most resistive part of the bar is clearly shrinking with depth, thus showing the 3D shape of the sand body. At larger depths (Slice 8—below 6.1 m) conductivity increases up to 180 mS/m, delimiting the base of the bar body. It must be noted, however, that the CMD Explorer provides, as a rule of thumb, reliable information only down to 6 m below ground and thus Slice 8 is effectively an extrapolation due to the need to have an infinite semispace at the bottom of the electrical conductivity model, and thus should be considered with care. Although the point-bar body shows a fairly homogeneous electrical resistivity, a subtle increase in resistivity values defines a 20 m narrow, NNE–SSW trending belt in the SE corner of Slices 2 to 6. The location of this belt fits with that of the abandoned channel forming the meander Bend 1 as apparent in satellite images (Figure 4b), and suggests that the higher resistivity values are linked to the coarser material of the deposits filling the abandoned channel. Deposits surrounding the low-resistivity point-bar body show conductivity values spanning from 80 to 250 mS/m, with values close to 100 mS/m down to 3 m below ground, increasing to 250 mS/m below 3 m. Comparison between geophysical data and geomorphic evidences suggest that these electrically conductive sediments represent floodplain deposits in which the Bend 1 meander was cut, thus developing the related point-bar sedimentary body.

3.3. Sedimentology

Core data provide the ground-truth information needed to calibrate/confirm geophysical data with localized information. The cores help define sedimentary features of the electrically resistive point-bar body, related channel-fill deposits and surrounding electrically conductive overbanks. All cores reveal that the deposits were intensely reworked by agricultural activities down to 80 cm below ground. Note that this fact may pose significant limitations to remote-sensing interpretation that is forcibly limited to surface images. Reworked deposits are dark brown and consist of very fine sand with a variable amount of mud. Point-bar deposits were completely cored at sites AV_a–c (Figure 7). Point-bar deposits occur from 0.8 m to a maximum of 5.4 m below ground and mainly consist of sand with a scarce percentage of mud. Cores AV_a–c reveal that bar deposits cover either organic-rich mud (core AV_b) or sandy deposits (cores AV_a and AV_c). Point-bar deposits are floored by a channel lag that consists of massive medium sand with pebble-sized mudclasts (Figure 8a). This basal lag is covered by lower bar deposits, consisting of 1–1.5 m of mud-free, well-sorted fine to medium sand, which is commonly massive or crudely plane-parallel stratified (Figure 8b). Upper bar deposits are ca. 2.5–3 m thick and consist of fine to very fine sand with subordinate mud layers. Sand is plane-parallel to ripple cross-laminated (Figure 8c) and contains mud for ca. 12%, 21% and 20% in the upstream, central and downstream zone, respectively. Mud layers (Figure 8c) range in thickness between 0.5 and 2 cm, and consist of massive or crudely laminated mud with plant debris. Lower bar deposits are ubiquitously mud free. The overall grain size of the bar deposits does not relevantly change along the bar, which appears as an almost monotonous sandy body from its upstream to downstream reach.

The overbank deposits were cored where geophysical investigations reveal the occurrence of electrically conductive sediments. These deposits mainly consist of silt-rich mud with subordinate sandy layers with horizontal bedding. Mud is massive and can be organic-rich (Figure 8d) or slightly oxidized in the lower and upper part of the overbank succession, respectively.

4. Discussion

4.1. Implications for Noninvasive (Remote Sensing and Ground-Based Electromagnetic) Investigations

This study shows how remote sensing and ground-based geophysical data represent an ideal combination of noninvasive techniques that can guide direct investigations and, at the same time, can be integrated to provide a 3D reconstruction of the shallow subsoil once supported by the local direct evidence for verification and calibration.

Our results show that in the study area the use of aerial images is very effective in supporting the rapid recognition of geomorphological and sedimentological features with a resolution approaching 0.5 m. However, the aerial pictures dating back to before the 1980s do not provide useful paleohydrographic evidence, probably because (a) conventional zenithal pictures were originally taken for cadastral purposes and, thus, were shot during the winter season when vegetation cover was limited; and (b) widespread leveling of the fields and use of strong plowing machines was common until that time, thus causing severe erosion of the topsoil, especially in the zones where convex landforms were formerly present, with slight accumulation in the depressions. Thus, in correspondence with natural levees and sand ridges related to scroll-bar sequences, the sandy and silty sediments were exhumed, showing a lighter signal in the soilmarks.

Cropmarks appear to be, in general, the most effective indicators of shallow geomorphological features in the considered environment. This is linked to the much higher permeability of the coarser sediments, leading to greater drainage and, in absence of irrigation, water stress to crops especially during the hot season (July and August). As documented in other areas of the Venetian Plain, e.g., [79], besides the weather condition of the period before the image acquisition, the maximum detail shown by the cropmarks strongly depends on the type of cultivated plants and, in particular, it decreases with the size of the leaves and the plant spacing. The best data are generally found in zones with soya bean or hay meadow crops. Wheat, barley and corn display variable visibility, as the first two are seeded in tight rows and have small leaves, but are harvested before mid July; in contrast, corn has a larger plant size and larger distance between rows, but is harvested in September or October and thus can (usefully) experience water stress. Note that this study was carried out analyzing freely available high-resolution images, reaching a resolution of about 0.5 m. This suggests that the use of specific images, acquired through latest commercial satellite or drone-borne multispectral scanners, could easily support the recognition of features between 0.1 and 0.5 m with superior results.

Geophysical data play a critical role in our analyses. They bridge the gap between surface-extensive information provided by remote sensing (that primarily guided the identification of the study area and the location of the surveys) with the information at depth carried out through traditional drilling and sampling operations. The latter in turn were positioned on the basis of remote sensing and geophysical evidence, thus minimizing the sampling to the locations where this information was needed for verification and calibration. The geophysical data constitute the backbone of the study in that they provide ultimately the 3D information to fully reconstruct the sedimentary structure of the site. This result, however, is not trivial to achieve.

First, great care must be posed in the acquisition phase. This entails not only the choice of the instrument and the strategy developed for covering a large area—in this case a (nonconductive) sledge towed by a tractor at sufficient distance not to induce current in the metal frame of the tractor itself. In addition, care must be paid when setting the instrument to have a good control of the measurement geometry, which is in turn essential to obtaining reliable inversion results. In this case, the stability of the sledge and the positioning care allowed us to obtain, for all data points, an RMS error of less than 10% between measured and simulated apparent conductivity data at the end of the inversion process. Note that carrying the same sonde by hand, particularly on rugged terrain, can easily unbalance the instrument, with measurements thus taken with some coils much closer to the ground than others. This might induce a very large measurement error that is impossible to correct a posteriori, as the true acquisition geometry is then completely unknown.

Second, obvious outliers must be removed from the dataset. These may include negative or extremely high apparent conductivity values (that are physically implausible for the given acquisition geometry), which may be due to local metallic or magnetic features.

Third, an accurate reconstruction of positioning must be made. The availability of reliable colocated data from HMD and VMD is essential to have the six independent pieces of information necessary for depth inversion. This requires both a guided pace in the field to reoccupy roughly the same locations and proper postprocessing to assemble the data that pertain to the same reasonable surrounding of each measurement point, in this case with a radius of 1 m. Noting that the sonde size itself is of a few meters (Figure 2a), this accuracy is perfectly acceptable for the purpose at hand.

Fourth, geophysical inversion is an inherently ill-posed problem. In this case, it means that while the general pattern of the electrical conductivity variation with depth is constrained by the data, the subtler details may not be retrieved uniquely. In other words, at each measurement location, a number of different 1D models, which all, however, maintain certain key patterns, can be equivalent in terms of goodness of fit to the measured data within the data uncertainty range. In particular, for example, if conductivity increases with depth, all inverted profiles will display the same general pattern. However, the transition from lower to higher conductivity may happen continuously, or in steps, and steps may occur at slightly different depths. While this can be viewed as a weakness of geophysical methods (and EMI in this particular case), it is not without remedy. Indeed, geophysics should never be applied without some “ground-truth” coming (as in this case) from drilling investigations (local and not without their uncertainties, but still necessary). Thus, an iterative revision of the inverted vertical profiles was conducted to select, among the plausible electrical conductivity profiles those that also fit reasonably to the direct evidence where this evidence is available. The procedure was simply performed manually, particularly selecting suitable layer interfaces compatible with drilling evidence. This type of procedure should always be applied, and it is in the energy exploration where 3D seismic data are blended, e.g., with well logs coming from deep borings. For the geophysical data used in this study, the “ground-truth” is given by the local comparison between lithology from sedimentary cores (AV_1, AV_a, AV_b and AV_c) and electrical conductivity derived from EMI inversion at the same locations. A plot showing the resulting correlation, also taking into account the variability of electrical conductivity for the same lithology (shown as standard deviation error bars) is shown in Figure 9.

4.2. Implications for Investigation of Meandering River Deposits

Integrated remote-sensing, geophysical and sedimentological approaches provide insights to discuss key features of investigated point-bar deposits, with specific emphasis on their genesis and internal grain-size variability.

The external boundary of this body is consistent with the morphology of a point bar that originated through lateral migration of a meandering channel [80,81,82]. Nevertheless, the curved profile of the inner boundary suggests that the bar accretion started from the inner bank of an arcuate channel that showed a sinuosity of ca. 2.3 (Figure 6b). This evidence contrasts the common assumption that point bars originate from a progressive increase in channel sinuosity of a relatively straight channel, which gradually migrates laterally until reaching a sinuous configuration [83,84,85,86]. The onset of accretion from a sinuous channel allows for storing a reduced volume of bar sediment in comparison with that produced by inception from a straight channel. In the study case, the 3D reconstruction of the bar body from EMI data shows that about 1.8 million m³ of sand are stored in the arcuate point bar body (Figure 10a). If one assumes that the accretion of this bar started from a straight channel, as would be suggested by the application of classical sedimentological models to remote-sensing data, the estimated volume of the bar would have been of ca. 3.1 million m³ of sand (Figure 10b), leading to a remarkable overestimation of the accumulated sand. The onset of point-bar bar accretion from a sinuous channel was probably driven by pre-existing floodplain morphologies, which, at the early stage of channel development, forced water to drain through the paths defined by adjacent depressed areas [87,88]. The establishment of a curved planiform morphology could have been forced by floodplain lithological and morphological heterogeneities [89,90], which are associated with the occurrence of numerous overbank subenvironments, including crevasse-splays, levees, floodplain lakes and floodbasins (cf., [91]). Different deposits and morphologies associated with these subenvironments possibly forced the newly formed watercourse to connect adjacent depressed areas and assume a sinuous geometry.

The location of the sedimentary cores within the point bar provides information concerning both vertical and lateral grain-size variability within this sedimentary body, with a particular focus on the comparison between upstream and downstream bar deposits. Although a fining-upward grain size distribution has been broadly considered to be typical of point-bar deposits [85,91,92], a certain variability of vertical grain-size distribution has also been documented [93,94]. Core data show that muddy layers occur in the upper part of the bar, although mud is visibly subordinate to sand. The grain size of sand varies significantly neither vertically nor downstream, and the bar is, therefore, characterized by widespread weak vertical grain-size trends. Although the lack of a clear vertical fining in the upstream bar zone is consistent with the occurrence of high bed shear stress [95,96], the paucity of muddy deposits in the downstream part of the bar is a peculiar feature, which cannot be ascribed to the overall lack of mud in the system, being that the overbank deposits are entirely made of mud. The open morphology of the bend [97] associated with the study bar could have hindered the formation of a dead zone, which commonly forms in sharp bends [18], preventing the accumulation of mud in the downstream bar zone.

Figure 11 shows a summary of the workflow we used to extract the relevant information from each data source and blend the pieces of information towards the final conceptual 3D distributed model of the studied site.

5. Conclusions

This paper presents a successful integrated approach to analyze the distribution of sedimentary facies of a paleomeander in the Southern Venetian Plain, Northeastern Italy. The approach is based on a combination of remote-sensing (aerial and satellite) data, geophysical investigations (electromagnetic surveys) and direct sedimentary coring.

From the methodological point of view, we show that the combined use of noninvasive techniques such as remote sensing and ground-based geophysical data provides an effective method for the purpose at hand. In particular, remote sensing is quite effective for the identification of sites of interest and features at a metric scale, which is potentially linked to different subsoil structures. In the case considered here, cropmarks are the most useful features observed from the satellite images, due specifically to the water stress induced in crops by the higher permeability of sandy bodies with respect to silty sediments. However, remote sensing can only provide information on the ground surface. On the contrary, geophysical methods are specifically designed to reconstruct the subsurface structure on the basis of contrasts of geophysical parameters. In this case, we used electromagnetic induction (EMI) methods, and particularly an FDEM small-scale multicoil system. Well-designed, acquired, processed, and inverted EMI data allowed us to extend the surface information provided by remote sensing to a maximum depth exceeding 6 m below ground level, allowing the construction of a 3D model of electrical conductivity of the subsoil. Direct investigations via sedimentary core drilling were positioned on the basis of remote sensing and geophysical data, in order to confirm and calibrate the geophysical investigations, which were also partly reinverted on the basis of the new evidence. The overall cycle of investigations thus allowed us to set up a 3D stratigraphic model of the site, consistent with all available data. On the other hand, the sequence of investigation activities was designed in such a manner that the information collected at one step optimized the design of the next step, thus reducing the overall effort required to complete the task.

From the sedimentary point of view, the point bar studied shows an uncommon arcuate morphology, that contrasts the common assumption that point bars originate from a progressive sinuosity increase of a relatively straight channel that migrates laterally until reaching a sinuous configuration. This can be explained by considering the variety of alluvial subenvironments in the floodplain. These floodplain heterogeneities likely controlled water fluxes over the platform, by facilitating water drainage within traces of depressed areas, defining the sinuous shape of the study channel during the very first phases of channel formation. As far as grain size distribution is concerned, although classical facies models highlight overall trends of upward and downstream fining of grain size within point-bar deposits, the grain-size trends of the study bar do not vary significantly either vertically or laterally. The bar is, indeed, characterized by a widespread weak vertical grain-size trend, and it appears as a homogeneous body of medium to fine sand. The lower mud content in the downstream portion was probably a result of the open morphology of the bend that could have prevented the formation of the dead zone, which is commonly directly linked to mud accumulation in the downstream portion of the sharp bends.

This study provides a solid basis for developing more detailed sedimentological investigations, which could be improved including acquisition of data concerning internal stratal architecture of the alluvial deposits. GPR investigations and recovery of undisturbed sedimentary cores would provide further relevant insights to this approach, with relevant follow up in the frame of subsurface exploration or management of surficial aquifers. Detection of the distinctive morphometric and sedimentological features of Late Holocene paleochannels would allow a comparison with those of the rivers draining the area currently, and allow quantification of human impact on riverine dynamics [18].