Petro-physical Characterization of the Shallow Sediments in a Coastal Area in NE Italy from the Integration of Active Seismic and Resistivity Data

Integrated geophysical analysis using different methods along with a priori information from wells, is a proven approach to investigate the geology and the petro-physical characteristics of subsoil. We collected seismic and geoelectric data in an area located on the Adriatic coast in North-Eastern Italy, aimed at characterizing the quaternary sediments and the shallow geological structures. Compressional and shear-wave data provided information about geometry and velocity of the quaternary sedimentary succession, while geoelectric data provided information about the resistivity in the shallower formation, which strongly depends on the presence of groundwater (brine) and on its salinity. Clustering analysis allowed us to subdivide the study area into subdomains showing similar values of resistivity and compressional- and shear-wave velocity, enabling for a better interpretation of the processed seismic sections. Then, we calculated the petro-physical properties of the investigated sediments, i.e., brine saturation and resistivity, porosity, and clay content, for each cluster. This inverse problem involves rock-physics theories and an optimization algorithm based on the simulated annealing global-search method. The results, validated using borehole stratigraphy, provided information about the salty water wedge intrusion.


Introduction
On the northern Adriatic coast, the balance between sea erosion and aggradation of the river sediments gives place to its iconic lagoons (e.g., Fontolan et al. 2007;Bondesan et al. 2004). Two main sources of sediments, on the beaches of the outermost islands which separate the lagoons from the open sea, can be identified. The first source, carrying the majority of the material, is the estuaries of the rivers, which provide relatively thickgrained sediments, which are then further transported by the coastal currents and form the characteristic long, sandy beaches. The second source, carrying finer-grained muds, rich in organic matter from the lagoons is the tidal inlets, which deposit sediments along fans at their edges (Escoffier 1940;O'Brien 1969;Fontolan et al. 2007).
Specifically, in this work we focus on a camping site beach close to the Bibione village, in the province of Venice, North-Eastern Italy (see Fig. 1). Previous geomorphological studies showed that the two main sources of sediments in the area are the Tagliamento river and the tidal inlet of the Porto Baseleghe lagoon (Bondesan and Meneghel 2004). Francese et al. (2014) also include an estimate of the stiffness of the terrain through surface-wave tomography. The authors deployed a set of seismometers in the whole Province of Venice and produced a map of the V S30 . In this map, the area appears to have a higher V S30 compared to the surrounding beach, probably due to higher clay content of the sediments compared to the eastern part of the beach, where the coarser sands of the Tagliamento river prevail.
In this study, we characterize the shallow subsurface to identify the sedimentary succession and the seawater-saturated medium with particular regard to the interface between fresh-and seawater. We propose a broad, integrated geophysical survey, in which both seismic and electrical resistivity data are collected and analyzed. Regarding the seismic data, we acquired and processed P-, SV-, and SH-data, as well as surface waves. Along the same line, we collected an electrical resistivity data. Such a comprehensive geophysical survey can rarely be found in the literature. Acquisition and processing of P-, SH-, and SV-wavefields have been mostly used for hydrocarbon exploration (Macrides and Kelamis 2000;Beecherl and Hardage 2004). Differences in SV-and SH-velocities are an indication of the presence of anisotropy, and the phenomenon of shear-wave splitting has been studied in a seismological context (Bale et al. 2009;Suroso et al. 2017). Applications of such concepts to the near-surface deposits are even rarer. An example could be found in Da Col et al. (2021), who performed a high-resolution seismic survey to characterize the flysch on the eastern Adriatic coast. In this work, we further integrate the results of highresolution seismic survey and of electrical resistivity tomography by performing a petrophysical inversion.
The term petro-physical inversion indicates an optimization procedure to obtain some physical properties of the rocks from the integration of different geophysical data (e.g., Sun and Li 2016;Kamm et al. 2015).
Integration of the different geophysical datasets using rock-physics theories can be performed either via a joint inversion (Sun and Li 2016;Lelièvre et al. 2012) or via a quantitative joint interpretation adopting different clustering methods (Wang et al. 2021;Bedrosian et al. 2007;Martinez and Li 2015). In both cases, the superior reliability and quality of the obtained results are well established in the literature, as data integration is the best way to obtain an unambiguous geological model, especially in geologically complex areas (Bosch 1999;Bosch et al. 2002;Colombo et al. 2008;Di Giuseppe et al. 2016). For example, Bedrosian et al. (2007) classify geological structures in the area of the Dead Sea from a joint parameter space including seismic and magnetotelluric data. Geophysical data integration can be very helpful in mining exploration. For instance, Martinez and Li (2015) apply a similar method to gravity and aeromagnetic data to identify a mineralized ore body in Brazil. Lelièvre et al. (2012) perform a joint inversion of seismic and gravity data by incorporating two petro-physical relations (one for the sediments and the other for the massive sulfide deposit) in the inversion scheme. Di Giuseppe et al. (2016) apply a k-means algorithm to resistivity, seismic, and density data to identify different lithologies and fluid saturations in a geologically very complex hydrothermal volcanic area in Southern Italy. Another application of geophysical data integration to a complex volcanic environment can be found in Garcìa-Yeguas et al. (2017). These authors apply a fuzzy c-means clustering method to magnetotelluric and seismic tomography data on the Canary Islands, which allowed them to identify a geothermal system about 600 m b.s.l., which is probably the source of the fumarolic activity visible on the volcano. Wang et al. (2021) present a workflow to perform clustering with a fuzzy c-means algorithm on frequency electromagnetic data, high-resolution magnetic, radiometric, and gravity data. Mollaret et al. (2020) propose a new petro-physical joint inversion scheme based on electrical resistivity and refraction seismic 2D profiles acquired in a permafrost area to estimate the near-surface water, ice, air, and rock distributions. Shahin et al. (2022) describe an integrated multiphysics technique to jointly model resistivity, density, and elastic responses of dual-pore sandstones from well-log data. They performed petro-physical inversion by employing a global optimization search engine to determine porosity, matrix properties, water salinity, and saturation.
A different approach consists in employing rock-physics templates (RPTs) to estimate the physical properties of the subsoil from the integration of different geophysical data (e.g., Picotti et al. 2018). For example, Pang et al. (2019) propose to use RPTs based on the 1 3 seismic velocities and attenuation (Q factor) to estimate porosity, saturation, and permeability from seismic data.
An analytical approach similar to that which we present in this work can be found in Benjumea et al. (2021) and Böhm et al. (2015). Benjumea et al. (2021) aim at detecting karstic voids below the quaternary sediments and characterizing their filling (clay content, fluid conductivities, and saturation). To achieve this purpose, they first invert the acquired surface wave and geoelectric data and then integrate the resulting shear-wave model with the electrical resistivity profile. A c-means fuzzy logic clustering algorithm is applied to detect four different lithologies. Then, to univocally detect the karstic voids, a P-wave velocity and resistivity dataspace is clustered. Böhm et al. (2015) performed an analogous petro-physical characterization from a synthetic cross-well seismic and EM dataset generated in a realistic aquifer, partially saturated with CO 2 . They used a forward optimization method to recover CO 2 saturation, porosity, and clay content using P-wave seismic velocity, seismic attenuation (Q factor), and electric conductivity.
In this work, we present a method which integrates seismic travel time and electrical resistivity data in an innovative way. First, we invert the first-break travel times and resistivity data to obtain P-and S-wave velocity and resistivity sections. Then, in order to identify media having similar characteristics both in lithology and fluid saturation, we perform a clustering on the dataspace containing the information from all three tomographic profiles. Finally, we apply a petro-physical inversion to each identified cluster, computing porosity, clay content, brine saturation, and resistivity. Such detailed characterizations of lose materials are generally performed in laboratory experiments, using ultrasounds and other testing machines (Zimmer 2003), or in situ instruments like piezometric bender elements (Yang et al. 2018) or core penetration tests (Lancellotta 2008). In this work, we obtain the sediments characteristics using a standard geophysical equipment, up to a depth of 30-40 m.

Geological Setting
The data were acquired on the western most part of the Bibione beach, very close to the Porto Baseleghe tidal inlet (Fig. 1). The main sources of sediments are the Tagliamento river, about 10 km to the east, and the tidal inlet of Porto Baseleghe. The volume of the tidal prism was strongly reduced by the reclamation of land for agricultural purposes; however, the prism volume remains considerable and equal to 3.3 × 10 6 m 3 (Fontolan et al. 2007). Bondesan et al. (2004) describe the granulometry of the sediments changing from sandy in the proximity of the river estuary, shifting gradually to clay and mud close to the tidal inlets. Generally, the size of the grains tends to decrease moving to the west. In fact, in the eastern part in selected areas gravel can be found or even, very rarely, pebbles as big as 1 cm. The Baseleghe inlet can be considered natural or almost natural, despite the dredging operations performed to guarantee the navigability. This, together with artificial deposits of sands, caused the westward aggradation of the sandy beach.
The camping site is protected from the sea swells and from the acqua alta by a natural dune barrier, approximately 2 m high. Furthermore, vegetation in the form of grass, bushes, and a few medium-sized trees, is present in the northern part of the line, starting from the sand dune.
The coring with sedimentological and paleontological analysis VV located at 2.7 km from the investigated area (see Fig. 1) evidenced the presence of post-LGM (Last Glacial Maximum) sediments (Olocene) down to a depth of 27.5 m. This sequence is characterized by engraving fill deposits between 10 and 27.5 m of an incised valley.
LGM sediments related to a continental environmental were found between 28 m and 31.5 (Upper Pleistocene) and are characterized by alternations of clayey, sandy, and peaty sills, with the presence of thin layers of fine sand. In addition, two erosive surfaces have also been identified at approximately 27 and 28 m depth.
Geological and geotechnical investigation in the study area evidenced that the water table depth is located within 1 m from the surface.
Even though boreholes are not available on the beach where the data were acquired, it is likely that, below the surface sand, clay muds from the lagoon, and part of the sediment fan coming from the tidal inlet are present. This is also confirmed by the V S30 computed by Francese et al. (2014), which show an increase in velocity from east to west along the beach. This confirms that clay sediments from the tidal inlet could be present below the recently deposited sands.
Finally, large amounts of algae are deposited on the southernmost part of the beach, where the swelling occurs during the winter storms. Every year a new layer of sand is deposited on top of these algae either by wind or artificially input to preserve the beach from sea erosion. The result is a layer, possibly a few meters thick, comprised an alternation of dry algae, and other organic material alternate with lose sand.

Seismic Data
The compressional, and longitudinal and transversal shear seismic data were acquired along the same line in June 2020 (P-waves) and in September 2020 (SH, SV, and surface waves), respectively. The profile is oriented north to south and is approximately 200 m long.
In the June 2020 survey, we adopted a receiver spacing of 4 m and a shot interval of 2 m. In the September 2020 survey, using more geophones, we adopted a receiver spacing of 2 m and a shot interval of 4 m. Recording time was 1 s with 1-ms sample rate for both surveys. Elastic energy was propagated in the subsurface using a vibrator, mounted on a wheelbarrow, able to vibrate both P-and S-waves. The S-wave vibrator was oriented perpendicular and parallel to the line direction to obtain SH-and SVwaves, respectively. Two sweeps were performed at each shot point and stacked.
To be noted that the source coupling, although not optimal for the presence of sand, was quite good thanks to the heavy rainfall which occurred the nights before the acquisition of both P-and S-waves.
The data were recorded using 10-Hz vertical geophones for the P-waves, and 4.5-Hz horizontal geophones for the S-waves. The horizontal geophones were oriented in-line for the SV acquisition and cross-line for the SH acquisition.
The surface-wave data were acquired using 48 vertical 4.5 Hz geophones, spaced 2 m, and using a weight drop as source. A couple of shots were collected at the beginning, at the center, and at the end of the line.

Geoelectric Data
The ERT profile was acquired in correspondence to the seismic line (see Fig. 1). The resistivity data were collected using a new concept multisource wireless data acquisition system (MPT-IRIS-www. iris-instr ument. com). This system (LaBreacque et al. 2013) comprises several stand-alone transceivers synchronized via GPS timing and controlled by a 900-MHz radio signal (wireless protocol). Each transceiver is capable of handling a maximum of three electrodes. is worth noticing that with this peculiar design, transmitting and receiving dipoles must belong to different units. The peculiar and innovative feature of this system, beyond the wireless communication between units, is its capability of transmitting the current simultaneously with multiple dipoles, each one connected to a different transceiver (Bocchia et al. 2021). These characteristics make this system suitable to work also in very difficult environments (Picotti et al. 2017).
The dataset consists of a total of approximately 3060 measures. The acquisition geometry consists in the dipole-dipole configuration, with the electrode spacing set to 5 m. The ERT profile was designed gathering the units in 3 blocks each one comprised 5 units ( Fig. 2), thus providing a maximum theoretical depth of investigation of about 40 m.
Particular attention was devoted to the coupling of the electrodes with the ground. Each electrode was watered using salt water to reduce the contact resistances and improve the data quality. Before every acquisition, a contact resistance check was done; after electrode watering, the measured values of the resistance were in the range of 200-800 Ohm.m for all adjacent electrode couples, indicating good operational conditions.

Seismic Data Processing
Two different processing procedures were applied to the acquired data. The first one was devoted to facilitate the picking of the first arrivals and therefore to maximize the reliability of the tomographic inversion, while the objective of the second one was to improve the signal/noise ratio and maximize the quality of the seismic imaging (Yilmaz 2001). The final time migrated section, obtained from the analysis of the reflected events, provides information about the geometry of the geological formation.
The processing sequence used for the picking of the first arrivals consisted in: • Application of a deconvolution before the cross-correlation of the data with the pilot traces (Baradello and Accaino 2016);

Fig. 2
Survey scheme adopted for the geoelectric acquisitions. Each block is comprised 5 units and each unit controls 3 electrodes. The block distance is computed using an electrode spacing of 5 m • Stacking of the shots made in the same positions, in order to increase the signal/noise ratio; • Application of a filter for the elimination of random and coherent noise; • Application of automatic gain control (AGC) to increase the signal amplitudes at large offsets.
The main source of random noise was human activities taking place close to the study area, even though very few tourists were present. A source of coherent noise at a constant frequency of 28 Hz was the pump for the water circulation in the swimming pool right next to the line. Wind and strong swells were additional sources of noise. Examples of processed P and SH shots are shown in Fig. 3.
The processing sequence used for the reflection seismic imaging was carried out with the aim to remove coherent and random noise. In particular on the pre-stack P-wave data were performed: • Zero-phase deconvolution; • Trimmed mean-dynamic dip filtering (TMDDF) to partially remove the ground roll (Holcombe and Wojslaw 1992); • Stacking velocity analysis on the pre-stack data sorted in CDP ensembles (Yilmaz 2001); • Stacking of the CDP gathers after the normal move out correction; • Time migration; • Application of a FX deconvolution to the migrated data to enhance the horizontal continuity of the reflected signal.
In the SH data, all these steps of processing were applied except the zero-phase deconvolution that did not provide adequate results.

First-Break Tomography
We performed a first-break tomography for all three wavefields (P, SH, SV) using the OGS CAT3D software (Böhm et al. 2014), which is based on a minimum time ray tracing (Böhm et al. 1999). This software uses the simultaneous iterative reconstruction technique (Stewart 1993). We picked the first breaks on the common-shot gathers as shown in Fig. 4. Picking at near offsets of the direct waves was not always reliable due to the interference with the ground roll, and in this case, as in the example in Fig. 4, it was not performed.
In the inversion procedure, we applied the staggered grid method (Vesnaver and Böhm 2000). This technique allows us to define the problem on a starting grid, with large pixels, where the null-space is low, and therefore the problem is well posed. Then, by shifting the grid several times in both X and Y directions, we obtain different grids and compute a separate inversion for each of them. We then sum and average the results of each inversion and obtain a much finer grid, still preserving the well-posed problem.
The final depth interval velocity sections obtained by the first-break tomography were then converted in time interval velocity by using the Seismic Un*x software (Stockwell 1999). The obtained results were then superimposed to the seismic sections to check the consistency between the position of the reflecting horizons and the velocity fields.

Surface-Wave Analysis
We performed the multichannel analysis of surface-waves (MASW) in three steps. First, we computed the f-k spectrum of the signal on the full array of 48 geophones (Fig. 5a). On such spectrum, we identified the most energetic, linear (i.e., dispersive) events, and identified the dispersion curve corresponding to the maxima within such event. We then removed from the curve points which were evident outliers. Figure 5b and c shows a plot of the final dispersion curve, respectively, in the frequency-phase velocity and frequency-wavelength domains. Finally, the dispersion curve was inverted to produce 1 million 1D shearwave velocity profiles using the neighborhood algorithm, provided by the Geopsy software (Wathelet 2008). The initial model uses three layers plus half space; the properties of each  Table 1. The values we use for density and Poisson's ratio are in accordance with the values for dry sand given in Zimmer (2003).
Based on the misfit values of the computed profiles, we then only considered the 1000 best fitting ones (as in, those generating the dispersion curve showing the lowest misfit with respect to the picked one) for interpretation.

Electrical Resistivity Tomography
The dataset was subjected to a quality control prior to the inversion procedure for the production of the resistivity model (imaging). Approximately 47% of the total quadrupoles were removed with the following criteria: removal of the reciprocal measurements; elimination of measurements with instrumental standard deviation (obtained through multiple measurement stacks) larger than 10%; removal of the receiver potentials (V) lower (in absolute value) than 0.01 mV; and filtering of negative apparent resistivity values. Figure 6 shows the histogram of the statistical distribution of potentials, currents (I), and apparent resistivities after filtering the noisy and the reciprocal measurements. The average injected current is 177 mA, while the average apparent resistivity is about 12 Ohm.m, which is a quite low value, reflecting the fact that the line is located close to the shoreline on both sides (see Fig. 1a).
The filtered data were then inverted using the ERTLab Studio software (http:// www. geost udias tier. it), in order to perform the imaging and characterize the correct resistivity distribution in the subsoil.
The three-dimensional inversion parameters were the following: Mesh size in the x, y, and z directions: 2.5 m, equal to half electrode spacing; Initial homogeneous resistivity equal to the average apparent resistivity value; Noise level of 1%, estimated by a statistical analysis of the reciprocal measurements. Figure 7a shows the progress of inversion where the vertical bars refer to the proceeding of Chi-square ( 2 ) statistic. The inversion converged after a maximum of 7 iterations, when the 2 equals the number of measurements of the whole dataset. Figure 7b shows the cross-plot of the modeled data versus measured data at the final iteration, evidencing the good agreement between modeled and recorded potential values.

Clustering Analysis
In order to identify sediments with different petro-physical properties, we perform a cluster analysis on the P velocity, SH velocity, and resistivity of the tomographic model, where the chosen parameters are grouped in classes (clusters) that contain homogeneous properties of each parameter.
Before performing the cluster analysis, we resampled the parameters in the same model with the same voxel discretization. We then considered only the common grid for the cluster computation. As for the cluster method, we used the k-means algorithm (MacQueen 1967), one of the most commonly used and reliable methods for cluster analysis. In this method, the number of classes must be determined "a priori". Thus, we tested several values and chose the number of classes that gave us the best distribution of zones with homogeneous characterization of parameters. We find this method to be more suitable for our purposes since, based on our experience, it allows a control of the results while changing the number of classes. The k-means algorithm is based on an iterative procedure that minimizes the variance of each group of parameters (classes) with respect to its centroid. The algorithm follows these steps: (i) extracting the data points of the parameters space corresponding to all the reliable grid points (pixels) of the model (only those interested by the inversions), (ii) normalization of all considered parameters to obtain a comparable spatial The inversion converged after 7 iterations, when the 2 equals the number of measurements of the whole dataset (horizontal red line). Cross-plot of the modeled data versus measured data at the final iteration (b), evidencing the good agreement between modeled and recorded potential values range of values for the clustering computation, (iii) choosing the number of classes. We tested several values and considering the best distribution of the zones with homogeneous parameters characterization. In this case, the tests indicate that 6 classes represent the best choice, (iv) preliminary estimation of the cluster centroids as a regular distribution in the parameter space, (v) assigning each data point to the nearest cluster centroid, (v) upgrading the centroid position as the new barycenter of the points associated with the same class, and (vi) repeating the two previous points until the new centroid positions do not change significantly.

Petro-Physical Analysis
Velocity and resistivity models derived from P-and S-wave seismic and ERT surveys can be integrated to infer the petro-physical properties of the subsoil. This inverse problem involves an optimization procedure to derive brine saturation S b and resistivity (inverse of conductivity) −1 b , porosity Φ , and clay content C . In this work, we follow a procedure similar to that of Böhm et al. (2015), with the difference that in this case we have the SH-wave velocity instead of the P-wave attenuation.

Forward Model
The theoretical approach used to model the seismic velocities is the generalized Gassmann's model for multi-mineral sediments, described in detail by Carcione et al. (2005). We consider a porous medium composed of two solids, sand (quartz) and clay, such that the dryrock bulk and shear moduli are given by K m = ∑ i K mi and m = ∑ i mi , i = 1, 2, where index 1 refers to sand (q) and index 2 refers to clay (c). A generalization of Krief's model for a multi-mineral porous medium is used to obtain the frame moduli, (Carcione et al. 2005(Carcione et al. , 2022, where K i and i are the solid grain bulk and shear moduli, K s and s are the Hashin-Shtrikmann averages (Hashin and Shtrikmann 1963;Carcione et al. 2022), i is the fraction of material i per unit solid volume, such that 1 = 1 − C and 2 = C . V K = ∑ i i K i and V = ∑ i i i are the bulk and shear Voigt averages, respectively. Moreover, A i and are dimensionless parameters, where A i depends on the pore shape and Poisson's ratio of the matrix, and is a fitting parameter. The parameter A i represents a pore compliance coefficient and takes a value of about 2 for spherical pores, increasing as the pores become more crack-like (Le Ravalec and Gueguen 1996;David and Zimmerman 2011;Carcione et al. 2005). The parameter takes into account other factors that might significantly affect the shear modulus, as the presence of organic matter in the soil (e.g., Ekwue 1990).
We obtain the effective bulk modulus of the gas (air)-liquid (brine) fluid mixture in the pore space using Wood's model (Mavko et al. 2009 where K a and K b are the air and brine bulk moduli. The effective sediment (bulk) density is simply given by the arithmetic average: = (1 − Φ) s + Φ f . The density of the grains is s = (1 − C) q + C c , where q and c are the sand-grain (quartz) and clay densities, respectively. The density of the fluid mixture in the pore space is f = 1 − S b a + S b b , where a and b are the air and brine density, respectively. Then, the P-and S-wave velocities are where the wet-rock shear modulus G is simply the modulus of the dry rock, m (second Gassmann's equation), and the wet-rock bulk modulus K G follows from Carcione et al. (2005Carcione et al. ( , 2022, where The complex refractive index method (CRIM) for a shaly sandstone with negligible permittivity and partially saturated with air, can be expressed as (e.g., Schön 1996;Carcione et al. 2012;Rossi et al. 2022), where q , c , b , and a are the sand grain (quartz), clay, brine, and gas (air) conductivities, respectively. Generally, gas conductivity can be neglected. If γ is a free parameter, the equation is termed Lichtnecker-Rother formula. It is based on the ray approximation. The travel time in each medium is inversely proportional to the electromagnetic velocity, which in turn is inversely proportional to the square root of the complex dielectric constant. At low frequencies, displacement currents can be neglected and the above equation is obtained. For zero clay content, and neglecting q and a , Eq. (5) is exactly Archie's law used in Hoversten et al. (2006). Figure 8 shows the seismic velocities V P , V S , and the resistivity −1 of a sediment, obtained by using Eqs. (2) and (5), respectively. The material properties of the sediment components are shown in Table 2, where the grain conductivities can be easily verified by using CRIM Eq. (5) and the conductivity measurements of fully water-saturated sediments (e.g., Knoll and Knight 1994). The curves in Fig. 8 are computed as a function of porosity Φ , by using C = 40%, S b = 99% (almost full brine saturation below the water table), b = 0.22 S/m (brackish water), and =1.5.

Global Optimization
The computation of the petro-physical parameters (namely Φ , C , −1 b , and eventually S b , if the sediments are partially saturated with brine and air), from the seismic and electrical geophysical properties, is based on the minimization of the following misfit function (e.g., Böhm et al. 2015), where W 1 and W 2 are weighting coefficients and V T P , V T S , and T are the tomographic P-and SH-wave velocities and conductivity, respectively. The minimization of function (6) is a nonlinear problem whose solution requires efficient derivative-free optimization algorithms based, for example, on the simulated annealing (Kirkpatrick et al. 1983) or pattern search methods (Griffin et al. 2008). In this work, a simulated annealing global-search technique was adopted, because it is particularly effective in optimization problems that involve functions with a large number of independent variables (Sen and Stoffa 2013). Simulated annealing is an optimization method able to distinguish between local and global minima of a multiparametric complex function that may have several optima. This method has been successfully used in a variety of geophysical studies (e.g., Denich et al.  2) and (5), respectively. The material properties of the sediment components are shown in Table 2, and the petro-physical properties are the following: C = 40%, S b = 99% (almost full brine saturation below the water table), b = 0.22 S/m (brackish water), and =1.5 Table 2 Typical material properties of a sediment composed by sand (quartz) and clay, whose pore space is partially saturated by brine and air 2023; Sen and Stoffa 2013, and references therein). The minimization process is based on the concept of temperature T, as a physical analogy of the cooling (annealing) of ideal crystals in thermodynamics (Kirkpatrick et al. 1983;Cerny 1985). After a cooling process, a crystal lattice reaches the equilibrium at a lower energy state. When minimizing a function, any downhill step is accepted, but also an uphill step may be accepted, in order to escape from a local minimum. During the annealing, each new solution is accepted according to a temperature-dependent probability function. The most common acceptance criterion is that of Metropolis (Metropolis et al. 1953). Annealing may start either from a randomly generated point or from a point decided by the user. At the beginning of the process, the higher T determines a larger probability of acceptance of new solutions, allowing the algorithm to escape from local minima and thereby to explore a wide region of the search space. As the optimization process proceeds, T reduces, the length of the steps decreases, and also the probability of acceptance of unfavorable solutions is reduced. The annealing procedure is iterated until convergence to the global minimum occurs.
In this work, we used the algorithm developed by Goffe et al. (1994) based on Corana et al. (1987). We successfully tested the method by using a variety of different starting models. Firstly, we computed the P-and SH-wave velocities of the sediments, as well as their electrical conductivity by using Eqs. (2) and (5), respectively. Then, we perturbed the analytical values by adding noise with a standard deviation of 1%. After the forward modeling, we performed the inversion by using the simulated annealing global-search technique (Goffe et al. 1994) with different weighting coefficients and optimization parameters. By computing the percentage deviations between the inverted and initial models, the tests indicate that the following parameters guarantee a good compromise between computational costs and inversion reliability: initial temperature of 10 6 , temperature reduction factor of 0.85, and error tolerance for termination of 10 −9 .
These performance tests are not only useful to set the parameters of the optimization algorithm but also for the choice of the weighting coefficient to balance the contributions of the different types of input data, i.e., seismic velocities and resistivity. Figure 9a displays the results of a test carried out using W 1 = W 2 =1 with the following initial model: Φ = 40%, C = 45%, S b = 95%, = 2.5 and different brine resistivity values from 0.2 to 20 Ohm.m. The sediment characteristics are shown in Table 2, while the adopted constraints are listed below: 1. Porosity ranges between 0 and 70%; 2. Brine resistivity ranges between 0.2 and 30 Ohm.m; 3. Clay content ranges between 0 and 70%; 4. Brine saturation ranges between 90 and 100% (almost full brine saturation below the water table); 5. The fitting parameter ranges between 2 and 3.
Then, to investigate the dependence of the parameter deviations on the weighting coefficients, we used different combinations of W 1 and W 2 . The tests show a significant decrease in deviations by increasing W 2 from 1 to an optimal value of 100 and no evidence of improvements by increasing W 1 . Higher values of W 2 do not lead to further significant improvements. Figure 9b displays the results obtained with the same initial model described above but using W 1 = 1 and W 2 = 100. Deviations (absolute values) decrease by an average factor of 1.6 from the former case of Fig. 9a to the latter case of Fig. 9b, while clay content and brine resistivity exhibit the largest errors in both cases. The maximum deviations for C and −1 b exceed 50% in Fig. 9a, while in Fig. 9b are lower than 20% and 48%, respectively. Furthermore, the tests allow to estimate the typical errors associated with the computed petro-physical parameters, for the type of sediments considered in this work. Adopting the above optimization parameters and using W 1 = 1 and W 2 = 100, the following average (absolute) errors have been estimated: 1.7% for Φ , 7.7% for C , 1.6% for S b , 15% for −1 b , and 1.5% for .

Results
The results of the tomographic inversion of the first arrivals of the P, SH, and SV waves are shown in Fig. 10. P-wave velocity field (Fig. 10a) highlighted the presence of a shallower velocity of about 500 m/s in the northern part and a velocity close to 1100 m/s in the southern part. Velocity of about 1300 m/s, between 3 and 9 m depth, is obtained in the northern part, while in the southern part, velocity values of 1300 m/s down to a depth of 14 m are detected. The tomographic inversion allowed to define, in the central part of the line, the velocity down to a depth of about 35 meters. In the deep part the velocity ranges between 1750 m/s and 2200 m/s.In order to check the reliability of the tomographic inversion a time residual analysis was performed, considering the difference between the picked travel times and the travel times computed from the final tomographic velocity model. For the final P velocity field an RMS value of 0.8 ms was obtained, corresponding to an error of 1.8% with respect to the picked travel times. SH-and SV-wave velocity fields (Fig. 10b and c) show a shallow layer with velocity and thickness of about 110 m/s and 1-2 m, respectively. Then, until a depth of about 27 m, the velocity increases to about 200-220 m/s. At depth greater than 27 m, the SV velocity Fig. 9 Performance test of the optimization algorithm using the following initial model: Φ = 40%, C = 45%, S b = 95%, and = 2.5. The material properties of the sediment components are displayed in Table 2. The test shows the deviations (errors) in the reconstruction of the initial model as a function of the brine resistivity −1 b , by perturbing the theoretical velocity and resistivity values of the sediments. The perturbation consists in adding noise with a standard deviation of 1%. The displayed errors are computed by using different weighting coefficients: W 1 = W 2 = 1 (a), and W 1 = 1 and W 2 = 100 (b). Clay content and brine resistivity exhibit the largest errors. Deviations reduce by an average factor of 1.6 from the former case (a) to the latter (b) values increase to 300-320 m/s, while at the same depth, the velocity field of the SH-waves highlights greater heterogeneity with velocities between 250 and 320 m/s. The calculated RMS of the time residual values for the SH velocity field was of 4.5 ms (error 1.73%), while for the SV velocity field the RMS value was of 4.8 ms (error 2.64%).
Finally, we perform the analysis of the ratio between the SV-and SH-wave velocities (see Fig. 10d) and the ratio between P-and SH-wave velocities (see Fig. 10e). The SV/ SH velocity section shows values around 1. Only in the shallowest southern part of the investigated area, anomalous values, around 1.6-1.8 were found.The P/SH velocity section shows greater heterogeneity at all depths. In fact, at shallower depths, values of 10-11 are present in the southern part and of 2-3 in the northern part, while at larger depths (10-20 m), an anomaly of minor entity is evident in the central part, which shows values of 7-8 versus values of 5-6 at the sides of the section.The time migrated sections obtained after the processing of the P and SH seismic data is shown in Figure 11. The tomographic velocity fields, converted from depth to time domain using Seismix Un*x software (Stockwell 1999), are superimposed to the sections in order verify the consistency between the results of the tomographic inversion and the seismic section.
The two seismic sections evidence a good agreement with the results of the first-break tomographic fields. The P-wave section shows a very shallow reflecting horizon in the northern part, while in the southern part no clear evidence of such interface is visible (label 1 in Figure 11a). A clear interface is present at about 50 ms (label 2 in Figure 11a). This reflector is clearly interrupted at a distance of 80 m with a step in the northern part of about 15 ms. This step agrees with the result of the tomographic velocity field.The SH-wave Fig. 10 Velocity sections from the first-break tomography for P, SH, and SV first arrivals. The SV/SH velocity and P/SH velocity are also shown velocity section evidences a very shallow wavy interface (label 1 in Figure 11b) and a clear reflection at about 300 ms (label 2 in Figure 11b). Also this reflector is abruptly interrupted at the distance of 80 m.
In Fig. 12, we show the 1000 best-fitting profiles computed by the neighborhood inversion algorithm. To be noted that the color of the line indicates the misfit, computed as in Wathelet (2008), the values of which are indicated in the color bar. Finally, the dashed red line indicates the maximum penetration depth and therefore the last velocity change which can be detected before the half space. The profile looks compatible with what is found by both the tomography and the reflection seismic profiles as the main features are detected. More precisely, the strongest similarities are those with the SV velocity profiles. This is not surprising, as the Rayleigh waves are the result of the interaction of P and SV waves with the free surface. More in detail, we identify the velocity changes at 3 m depth, 12 m depth and, the sharpest of all, at 27 m depth. This is the last interface detectable with this method, since the maximum available wavelength is 57 m, therefore the maximum penetration depth (λ/2) is 28 m.
The result of the inversion of the geoelectric data is displayed in Fig. 13. As expected, the ERT model shows a general decreasing trend versus depth of the resistivity values, because the line is located close to the shoreline. Accordingly with the Ghyben-Herzberg relationship (e.g., De Wiest 1965), the fresh/salt water interface depth should be approximately 40 times the difference between the mean sea level and the fresh water table elevation. Therefore, the low average resistivity value at about 13 m depth likely indicates the presence of salt water, while the shallow high resistivity zones probably indicate the presence of three areas with dry sand deposits at near surface.

Discussion
Analysis of the results allows us to obtain new information about the geometry of the main geological horizons of the quaternary sediments, as well as their geological and petrophysical properties.

Geology
The investigated area is characterized by the presence of the megafan of the Tagliamento river. In this area, several transgressive incised valleys developed, and about 5 km of prograding sand ridges is observed along the coastline (Amorosi et al. 2008). In Fig. 14, a geological section until a depth of about 20-25 m, obtained from cores acquired in the area, is shown (see Fig. 1b for the position of the section). This section evidenced the presence of an incised valley, filled of sediments in correspondence of Concordia. Moreover, close to Caorle, it shows a coastal wedge overlapping the LGM alluvial plain that was dated, through radiocarbon method to 6 to 6.5 ka BP (Amorosi et. al 2008).
The base of the LGM alluvial sediments was defined by the analysis of a core acquired close to the investigated area (see Fig. 1b for the position of the core, VV). This core evidenced that the base of LGM is at about 31 m below the surface (see Fig. 15).
Multi-method geophysical imaging provided several information on the properties of the superficial sediments and on the geology down to a depth of about 30-40 m. Furthermore, the agreement between the outcome from different methodologies allows us to be confident about the reliability of results.
The comparison between the P-and SH-wave tomography and the corresponding seismic sections (see Fig. 11) shows a good consistency. In particular, concerning the P-waves, the interface placed at about 0.05 s (about 10-12 m) is well highlighted with both the reflection seismic and the tomography.
From the geological point of view, the discontinuity at about 0.05 s is the basis of the post-LGM deposits  and it is an important geological marker of the NE Italian plain, representing the passage between the Pleistocene and the Holocene. The depth estimated by our surveys is in agreement with the depth in the Caorle geological section (see Fig. 14).
As previously mentioned, in areas like the one considered in this work (post-LGM in the Veneto-Friuli plain), where watercourses flowed into the Holocene sediments, the situation is rather complex, as several transgressive-regressive cycles have occurred (Fontana 2006). An example of these cycles could be represented in the resistivity section where it is possible to identify three surface areas with high resistivity. These areas are probably related to sand dunes formed during different periods of marine transgression-regression.
Another significant result of the glaciers melting is the erosion of Pleistocene sediments in areas with flow of water. An example of this phenomenon can be found along the line  This erosion is evident both in the seismic section and in the velocity field obtained from the inversion of the P-wave first arrivals. The presence of these incisions agrees, for example, with Fontana et al. (2010a, b).
The superposition of the SH seismic section with the first-break tomographic velocity field of the SH-waves (see Fig. 11) highlights the presence of an important interface at about 0.3 s (corresponding to 27 m). This discontinuity is also evident in the results obtained from the inversion of the surface waves (see Fig. 12), allowing to be confident about the reliability of these results. This interface, in accordance with the results of the analyzed cores of the well VV (see Fig. 1b) for the position, can be interpreted as the base of the LGM sedimentary sequence or as an erosional surface close to the base of the LGM sedimentary sequence.
The difference in the depth of the top of the LGM deposits between well VV and the section in Fig. 14 is due to the fact that the VV well is located on an incision of the LGM deposits and then filled with Pleistocene deposits. The top of the incision in the VV well is placed at 10 m in depth (Fontana et al. 2010a), in accordance with the depth of the roof in the seismic line identified in section P.
The analysis of the SV/SH velocity section (Fig. 10d) and P/SH velocity section (Fig. 10e) evidenced in the southernmost part of the profiles, i.e., that closest to the Fig. 13 Results of the electrical resistivity tomography. The inverted resistivity section shows a general decreasing trend of the resistivity values versus depth. Accordingly with the Ghyben-Herzberg relationship (e.g., De Wiest 1965), the low average resistivity value at about 13 m depth likely indicates the presence of salt water seashore, that the shallow sediments (up to ~ 5 m depth) present the higher values of the sections. This is most likely due to the presence of alternating layers of algae and other organic material with layers of sand. In fact, in such area, the algae are deposited during the winter storms and dried by the sun. These are then covered by a layer of sand both naturally and artificially (with the summer season approaching, the tourist operators cover such algae with sand). Such layering, together with the intrinsic anisotropic behavior of organic material (the closest examples in the literature are related to peat soils (e.g., Landva and Pheeney 1980;Hendry et al. 2012)) gives place to strong anisotropy. Similarly, the strong difference in bulk and shear moduli of such materials, together with its low consolidation, leads to the extremely high Vp/Vs. Further investigations along with laboratory tests are mandatory to assess this anomalous soil behavior.

Clustering
In Figure 16(a, b, and c), we show the 3D display of the parameters distribution and the associated classes obtained from the cluster computation. The obtained classes were then assigned to the corresponding pixels of the model (Fig. 16d) in order to identify the areas of the model characterized by similar properties.
In this case, class 6 (high P and SH velocities associated with low resistivity) is located on the lower part of the model; anomalous high resistivity (class 5) is present in a small shallower part, north of the line; meanwhile, low SH velocities with high P/SH velocity values (class 2) appear on the shallower southern part of the line. It is important to underline that the two geological markers, top and bottom of the LGM alluvial sediments, are well identified and can distinguish between sediments belonging to different investigated formations, showing therefore different petro-physical properties.

Petro-Physical Properties
We adopted the procedure described in Sect. 4.6 to compute the petro-physical properties associated with the classes identified by the clustering analysis. To this aim, it is necessary to model the dry-rock and wet-rock moduli of the sediments, as well as their electrical conductivity. We consider a porous medium composed of two solids, sand (quartz) and clay, partially saturated by brine and air, whose typical material properties are shown in Table 2. The conductivity is computed with CRIM Eq. (5), while the dry-rock and wet-rock moduli follow from Eqs. (1) and from the generalized Gassmann's theory (Carcione et al. 2005(Carcione et al. , 2022, respectively. The computation of the petro-physical properties of the sediments (i.e., brine saturation S b and resistivity −1 b , porosity Φ , and clay content C ) requires an optimization technique to directly solve the minimization problem consisting in finding the global minimum of misfit function (6). Due to the complexity of this function, a derivativefree minimization algorithm based on the simulated annealing global-search technique (Kirkpatrick et al. 1983;Goffe et al. 1994) was adopted. The optimization algorithm, described in Sect. 4.6, was applied to the centroid of each clustering class, using an initial temperature of 10 6 , a temperature reduction factor of 0.85, and an error tolerance for termination of 10 −9 . In this case, as explained in Sect. 4.6.2, the temperature is an adimensional parameter. Moreover, in misfit function (6), we used the following weighting coefficients: W 1 = 1 and W 2 = 100. As described in Sect. 4.6, these optimization parameters and weights guarantee a good compromise between computational costs and inversion reliability (e.g., Goffe et al. 1994).
The following constraints were adopted: 1. Porosity ranges between 10 and 70%. 2. Brine resistivity ranges between 0.2 and 60 Ohm.m. 3. Clay content ranges between 0 and 100%. The maximum clay content is set to 20% close to the surface, because direct observations indicate that the topsoil sediments are mainly composed by sand. 4. Brine saturation ranges between 90 and 100% below the water table, located within about 1 m depth from the surface, while it ranges between 0 and 90% above the water table. 5. The parameter must be properly constrained, in order to allow the algorithm to correctly fit both seismic velocities V P and V S . As explained in the point 6 below, this parameter can be considered as an indicator of organic matter content in the soil. It ranges between 1 and 2 for class 1 and between 2 and 3 for all the other classes.
The results, shown in Table 3, are also displayed in Fig. 17. These results are in agreement with the local borehole stratigraphy (located along the coast, between Bibione and the Tagliamento river; ISPRA 2021) and with geophysical and geochemical investigations on the salt water contamination in the adjacent Venice Lagoon mainland (Di Sipio et al. 2006).
We highlight the following points: 1. Porosity decreases significantly with depth, due to compaction. 2. Class 1 shows higher clay content and brine resistivity than class 2. This result is consistent with the fact that, being class 2 located very close to the shoreline, its sand component is predominant. Moreover, brine salinity increases and −1 b decreases approaching to the sea. The higher porosity of class 1 compared to class 2 is probably due to the difference in (unconsolidated) clay content. Table 3 Input seismic-electrical parameters to the global optimization algorithm (left) and results of the petro-physical inversion (right). The input data are the centroids of the classes identified by the clustering analysis. The sediments and brine resistivities are −1 and −1 b , respectively. The different classes are ordered from the shallowest (class 1) to the deepest (class 6)

Seismic-electrical tomographic parameters
Petro-physical parameters  (Palacky 1987). Therefore, for the deepest classes 3 and 6, the sediments are fully saturated by low-resistivity brackish water. Di Sipio et al. (2006) found similar results along the Adriatic coast between Bibione and Venice. 4. Class 3 shows higher clay content than all the other classes, in agreement with the local stratigraphies (ISPRA 2021). Because of compaction, class 3 has higher porosity than class 6. 5. Class 5 has somewhat anomalous values: the clay content is almost negligible, the brine saturation is about 60%, and the brine resistivity is about 58 Ohm.m. These values suggest the presence of a sand dune which, being close to the surface, it is partially saturated with high-resistivity fresh water and air. This anomalous resistivity value might perhaps also be associated with the presence of algae. 6. Variation of the parameter between different classes probably depends on the presence of organic matter in the sediments. In particular, occurrence of peat reduces the shear strength of the soil (i.e., Ekwue 1990). The study area is characterized by algae deposits along the shoreline and by the presence of peat in the subsoil, as reported in the local stratigraphies (ISPRA 2021). This evidence might explain why in our case > 1.
Petro-physical results highlight a large difference in clay content of class 3 with respect to the other classes (see Table 3 and Fig. 17b). This class is attributable to sedimentary deposits of the LGM phase. In this phase, as a consequence of the glacial advance, the sea level dropped and the plain extended to the south for 400 km more than the present, showing a very low slope angle (less than 1‰). This low slope meant that the fluvio-glacial currents were not confined laterally in incision valleys, losing the ability to drag sands after 20-25 km from the glacial front (Fontana et al. 2010a). Only silts and clayey silts were deposited in the areas where the coast is currently present. This phenomenon can well explain the high percentage of clay obtained by the petro-physical inversion for class 3 with respect to the other classes attributed to interglacial sediments.
Another important result obtained by petro-physical inversion is the presence of salt water in the investigated formations below the sea floor (see Fig. 17c). The phenomenon Fig. 17 Results of the petro-physical analysis described in Table 3 and related to the classes identified by the clustering analysis displayed in Fig. 16 of saline intrusion is well known along the north-eastern Adriatic Sea. There are studies in the south of Venice (Carbognin et al. 2005;Di Sipio et al. 2006) which highlight how the saline intrusion can extend to several km inland reaching even great depths. This phenomenon, as far as surface structures are concerned, is affected by the morphology and by the presence of canals or rivers along which salt water can rise again, for example, in times of drought.
A report on the salinity of soils (Municipality of Caorle, 2014) is available for the study area. This report highlights the increase in salinity starting from 100 cm below the ground, in accordance with the resistivity section, and with the results obtained by the petro-physical inversion.

Conclusions
The integrated analysis of results from different geophysical methods, supported by a priori well information, allowed us to characterize the quaternary sediments of the coastal area close to Bibione (North-Eastern Italy), both from the geological and petrophysical point of view. In particular, analysis of compressional and shear seismic waves provided information about the geometry and the depths of the main geological markers, such as the passage from Holocene to Pleistocene and the bottom of the Last Glacial Maximum (LGM).
Tomographic inversion of first arrivals and electrical resistivity tomography gave information about the compressional-and shear-wave velocities, as well as the resistivity of the sediments down to a depth of about 35 m. Clustering analysis enabled us to define 6 classes of velocity and resistivity values along the investigated profile, yielding a useful integrated model for the geophysical and geological interpretation of the subsurface stratigraphy. The surface-wave data inversion validated the results obtained from the first-break tomography of the SH seismic data, confirming the depth of the LGM sediment base (about 30 m). This result is also in agreement with the stratigraphy of a local borehole.
The analysis of the resistivity section provided information about the shoreline location during different periods of marine transgression-regression events. These locations are characterized by the presence of sand dunes showing high resistivity values.
Finally, the implementation of an optimization algorithm based on rock-physics theories and on the simulated annealing global-search method allowed to further integrate the different data and to calculate the petro-physical properties of the investigated sediments, i.e., brine saturation and resistivity, porosity, and clay content, for each clustering class. The results are in agreement with local borehole stratigraphy and confirmed the high clay content characterizing the sedimentary deposits of the LGM phase. Moreover, the brine resistivity decreases versus depth, evidencing that the salt-wedge intrusion, at about 100 m inland, is located at about 20 m depth.
These results are important in the recent climate change context because coastal areas could be subjected to significant impacts due to changes in sea level rise, rainfall regimes, and temperatures, which in turn modify the salt-wedge intrusion, altering soil and inland water characteristics.
The petro-physical procedure described in this work could be successfully applied to other study cases, as for example, the evaluation of the integrity of river levees and earth dams.