Hydrogeology Journal (2014) 22: 691–703 DOI 10.1007/s10040-013-1087-8 A three-dimensional hydrogeological–geophysical model of a multi-layered aquifer in the coastal alluvial plain of Sarno River (southern Italy) R. Di Maio & S. Fabbrocino & G. Forte & E. Piegari Abstract The coastal alluvial plain of Sarno River (Campania Region, southern Italy) is a very rich environment that has experienced a long history of changes due to both natural phenomena such as eustatic sea-level variations and deposition of volcanoclastic sediments, and human civilizations who populated this area since historical times. As a result, it is characterized by complex stratigraphic sequences and groundwater flow systems. The architecture of the multilayered aquifer system in a sample area, located in a densely urbanized sector at the mouth of Sarno River, was reconstructed. Starting from the analysis of stratigraphic log data and laboratory geotechnical measurements, the lithostratigraphical-unit sequence was retrieved and a realistic three-dimensional (3D) model of the hydrogeological heterogeneity was obtained. The results of a detailed 2D electrical resistivity tomography survey were used to support the analysis of the spatial heterogeneity of the aquifer system in a sector characterized by lack of log data. The integration of hydrogeological and geophysical data allowed for the reconstruction of a 3D hydrogeophysical model of the multilayered system, which electrically characterizes and geometrically identifies two aquifers. Finally, piezometric-level Received: 15 May 2013 / Accepted: 22 November 2013 Published online: 19 December 2013 * Springer-Verlag Berlin Heidelberg 2013 R. Di Maio : S. Fabbrocino : E. Piegari ()) Department of Earth Sciences, Environment and Resources, University of Naples Federico II, Naples, Italy e-mail: esterpiegari@gmail.com R. Di Maio e-mail: rodimaio@unina.it S. Fabbrocino e-mail: silvia.fabbrocino@unina.it G. Forte Structural and Geotechnical Dynamics Laboratory StreGa, University of Molise, Termoli, Italy G. Forte e-mail: giovanni.forte@unina.it measurements validated the hydrogeological–geophysical model and showed the effectiveness of the methodology. Keywords Heterogeneity . Groundwater management . Electrical resistivity tomography . Hydrogeophysical modeling . Italy Introduction The full assessment of complex, heterogeneous and anisotropic aquifer systems is key to understanding groundwater conditions and contaminant transport. Groundwater models are useful tools for management and protection of water resources, but their predictive function is often limited by a poor characterization of hydrogeological parameters at a local scale. A variety of structural, volcanic and depositional processes produces a composite spatial distribution of hydraulic conductivity in alluvial aquifer systems. Groundwater flow and solute transport simulations in similar heterogeneous bodies of intercalated aquifers and aquitards need a large amount of observed data (Felletti et al. 2006; Ouellon et al. 2008; Hsien-Tsung et al. 2010; Vienken and Dietrich 2011). Note that in the following, for practical purposes, aquitards and aquicludes among confining units are distinguished according to Poland et al. (1972). Specifically, an aquiclude is considered as a body of saturated but relatively impermeable material, which has very low values of “leakance” (i.e. the ratio of vertical conductivity to thickness) and allows negligible interaquifer flow. An aquitard is a saturated poorly permeable bed, whose leakance ranges from relatively low to relatively high values (Laney and Davidson 1986). Many tools to acquire the knowledge of lithologic and stratigraphic features and spatial heterogeneity at different scales are available (de Marsily et al. 2005; Galloway 2010), but it is essential to improve the hydrofacies models including in situ geophysical investigations that better handle the rock strata connectivity and the geometry of aquifers. In such a context, electrical geophysical exploration may provide a reliable model of aquifer heterogeneity, as it is able to identify the strata sequence in terms of geometrical and physical parameters (Drahor et al. 2011 and references therein). In addition, proper relationships between aquifer electrical and hydraulic properties can be recognized in order 692 to define the groundwater potential (Kelly and Frohlich 1985; Huntley 1986; Ahmed et al. 1988; Chandra et al. 2008). Nevertheless, established standard methodologies for the assessment of hydrodynamic parameters of complex aquifer systems are not yet available (Linde et al. 2006; Mele et al. 2012). During the last few decades, different approaches have been proposed to correlate hydrogeological and geophysical data (Hinnel et al. 2010 and reference therein). These approaches can be mainly divided into three categories. The first two categories integrate geophysical and hydrogeological measurements by means of an uncoupled inverse approach, which relies on an independent inversion of geophysical and hydrogeological data. In this framework, it is possible to distinguish procedures that infer the hydrogeologic model by (1) simple inversion (e.g. Cassiani et al. 1998; Hubbard et al. 1999; Tronicke and Holliger 2005; Koch et al. 2009; Sinha et al. 2009) or (2) joint inversion (e.g. Hyndman et al. 1994; Chen et al. 2006; Linde et al. 2006) of survey data. In both cases, although different methods are proposed to reduce the uncertainties related to the non-uniqueness of the inversion procedure, geophysical and hydrogeological data are treated independently, and the hydrogeological model is not used to constrain the geophysical interpretation. Recently, a third category of approaches is proposed which relies on direct coupling of hydrologic models and geophysical models during inversion (Kowalsky et al. 2005; Ferré et al. 2009; Hinnel et al. 2010). In this study, an uncoupled inverse approach is adopted to obtain two three-dimensional (3D) models, respectively, hydrogeologic and geoelectric, derived from data sets acquired on two adjacent areas. Specifically, the geoelectrical survey was carried out to expand the hydrogeological model in a restricted area, where stratigraphic log measurements were not allowed. Then, a 3D hydrogeophysical model is set up by using the electrical resistivity inverted data as constrains for the analysis of the spatial heterogeneity of the aquifer system. The test area is the coastal alluvial aquifer system of the Sarno River plain (southern Italy). Despite its extension and complex geological history, the current stratigraphical and hydrogeological information is heterogeneous both for quality and spatial distribution of data. Due to the composite hydrodynamic condition and the high human impact of this plain area, the study is focused on the portion including the Sarno River mouth. In the following, after a brief description of the geological and hydrogeological setting of the selected area, its hydrostratigraphic and geophysical models are discussed and combined to characterize at field scale the multi-layer aquifer. Geological and hydrogeological features of the test area The test area is located at the mouth of the Sarno River plain (Campania Region, southern Italy), which is a coastal plain about 200 km2 in area and is bounded by Mesozoic carbonate reliefs to the east and south, and the SommaVesuvius volcano to the north (Fig. 1). It is the southern part Hydrogeology Journal (2014) 22: 691–703 of the peri-Tyrrhenian Campanian Plain Graben, originated during the Plio-Pleistocene by the opening of the Tyrrhenian Sea (Patacca et al. 1990). Its buried structure is composed of step-faults blocks of Mesozoic carbonate units of the Apennine chain (Ortolani and Aprile 1985; Brancaccio et al. 1991). NW–SE and NE–SW normal faults border the carbonate bedrock that dips towards the Somma-Vesuvio volcano and reaches the maximum depth of 2 km below the town of Pompei (Barberi et al. 1980; Di Maio et al. 1998; Gasparini 1998; Zollo et al. 1996, 2001), though it outcrops seaward as the limestone horst of Rovigliano islet, just close to the mouth of the Sarno River (La Torre et al. 1983). According to the geological evolution of the Campanian plain, the Sarno River plain experienced periods of tectonic subsidence and ground movements due to the volcanic activity and eustatic sea-level changes during the Quaternary (Brancaccio et al. 1991; Pescatore et al. 2001; Rolandi et al. 2003; Marturano et al. 2011). As a result, in this area, the Quaternary sedimentary successions are more than 2,000 m thick and their stratigraphic framework is very complex even at basin scale. Many studies, based on stratigraphic and micropalaeontological analyses of boreholes drilled in the south-western sector of the Sarno River coastal plain (Pescatore et al. 2001; Marturano et al. 2011), allowed the reconstruction of its paleoenvironmental conditions and sedimentary facies changes related to volcanoclastic aggradation/progradation and seaward shift of the shoreline. At basin scale, tephra deposits of the Campanian Volcanic Zone (Rolandi et al. 2003), namely of the Somma-Vesuvius volcano in the last 25 ka (Santacroce et al. 2008 and references therein), interbedded with marine and continental ones, provide important stratigraphical and geochronological markers. The upper part of the fill overlies an ignimbrite formation (Campanian Ignimbrite, 39 ka) related to one of the largest eruptions that occurred in the Mediterranean region in the last 200 ka, which vented along regional faults and uniformly covered the plain (Barberi et al. 1978; De Vivo et al. 2001; Rolandi et al. 2003). The tectonic pattern and erosion phenomena affect the toplap surface of its tufaceous facies, the so-called Campanian Grey Tuff, at depths of 10 to more than 30 m below sea level (b.s.l.). So that, because of tectonic dislocation and erosion phenomena, it is missing just at the mouth of Sarno River (Aprile et al. 2004). In summary, Pliocene and Quaternary deposits of the Sarno basin are both marine and continental as complex interactions between depositional processes evolve spatially and temporally. As demonstrated by deep boreholes, such as the Trecase 1, 1,900 m deep (Balducci et al. 1983; Brocchini et al. 2001), the Sarno Plain was above sea level during the Pliocene (Brancaccio et al. 1991; Cinque et al. 1993). From the beginning of the lower Pleistocene, subsidence phenomena occurred and caused marine ingression. In the lower part of the late Pleistocene, the subsidence, coupled to the volcanic activity, continued. Then, in the upper part, the plain uplifted due to the subsidence rate decrease and due to the last glacial regression. After the Campanian ignimbrite eruption and the consequent DOI 10.1007/s10040-013-1087-8 693 Fig. 1 Geological and structural sketch of the Sarno River plain (Campania Region, southern Italy) volcanoclastic aggradation, the plain became tectonically rather steady and fluvial down-cutting also took place for the sea-level lowering of the last glacial regression (Brancaccio et al. 1991; Romano et al. 1994). In the Holocene transgression maximum, lagoon and swamp systems have been detected more than 2 km inland from the present shoreline (Barra et al. 1996). As a result, significant spatial variations of permeability are observed. The groundwater flow occurs in the pattern of permeable coarse-grained sediments. Accordingly, more overlying groundwater flows occur, though at basin scale only one groundwater flow system was detected due to natural hydrogeologic conditions and well completion (Celico and Piscopo 1995; Fabbrocino et al. 2007). However, leakage phenomena among overlying aquifers can take place and the discharge occurs into the sea:first the river is gaining, while towards the mouth it is losing. 3D Hydrofacies model 3D geological and hydrogeological modelling is a composite process of learning, interpreting and visualizing buried depositional architectures, which is constrained by limitations of actual understanding and data availability (Robins et al. 2005; Beven 2007). A reliable modelling methodology has to recognize the value of different types of data in terms of both time and space as well as to make the existing data (including past model predictions) Hydrogeology Journal (2014) 22: 691–703 available. In this regard, it is possible to learn about places and to constrain predictive uncertainties by extrapolating information in sites where the data are not available and investigation cannot be allowed (Beven 2007). Following such criteria, this report provides the reconstruction of the hydrogeological block diagram for a complex coastal alluvial plain by combining the most common existing data for geological and hydrogeological characterization with high-resolution 3D resistivity modelling. The sample area extends for 250,000 m2 at the Sarno River mouth and its digital model derives from fit survey data. The study focuses on the integration of existing information, which consists of 14 stratigraphic logs at depths of 3 to more than 20 m, mainly located in the SE sector (Fig. 2), and some geotechnical parameters of recognized alluvial, marine and pyroclastic deposits derived from laboratory tests. The reliable reconstruction of the aquifer model is carried out at the scale of aquifer systems and hydrogeological complexes. The hydrostratigraphy characterization derives from the analysis of the genetic stratigraphy rank of depositional systems and system tracts (Catuneanu 2006), and involves hydrogeological survey and examination of cartographic records (aerial photo, geologic and topographic maps). Taking into account the Holocene evolution of the Sarno River coastal plain (see section ‘Geological and hydrogeological features of the test area’), the interpretation of stratigraphic logs and the identification of the depositional DOI 10.1007/s10040-013-1087-8 694 Fig. 2 Geological map of the sample area with distribution of existing stratigraphic logs facies arise from the description of the lithology as well as the grain-size composition and sorting of sediments. At reference scale, the sedimentological analysis of 14 available stratigraphic logs was focused on detectable changes in lithology, color and texture. Then, grain size, shape and sorting analyses were performed to identify geologic units related to specific volcanic events and progradation of the plain (Pescatore et al. 2001). Such a reconstruction of lithofacies evolutive trends highlighted relevant properties of sediments, resulting in permeability changes as well as geometric elements constraining hydrofacies. Figure 3, for instance, shows the hydrofacies identification process based on the aforementioned stratigraphic correlation for the stratigraphic-log B5 (see Fig. 2). Table 1 reports the hydraulic conductivity, K, obtained by falling head permeability tests of nine hydrostratigraphic units detected in the investigated subsurface. Borehole data were organized into segments and contacts. Each contact has been associated with a particular horizon, and therefore to the depositional sequence (Lemon and Jones 2003). The 3D hydrofacies model of the examined aquifer system was obtained by coupling hydrostratigraphic cross-sections with the main geological constraints correlated to geomorphological history of the plain. Figure 4 shows the hydrogeological fence diagram that is the framework of the 3D hydrofacies model obtained by using the Aquaveo GMS 8.0 software. The close shoreline suggests the site should be divided into two domains traced by the coastal onlap. Thus, beach and marine deposits are close to the current coast line, Hydrogeology Journal (2014) 22: 691–703 while alluvial sequences are in the inner plain as a result of interbedding and lateral facies changes (see Figs. 2 and 4). The beach lithofacies is mainly made of well-sorted dark yellowish and dark greyish sands from coarse to fine, with rounded clasts of lapilli and fragments of shells. It is representative of a typical shoreface facies and constitutes the Marine Sands Complex (MSC), which is prevalent in the western sector. The alluvial sequence constitutes a composite aquifer system (see Table 1 and Fig. 4), whose borehole-derived lithofacies were related to the widespread units (see Fig. 3) reported by Pescatore et al. (2001). The deeper hydrogeological complex, i.e. the Sandy Gravelly Complex (SGC), dates back to pre-Roman and Roman time, and is characterized by a lateral facies change from shoreface to foreshore facies, with grain-size coarsening. The SGC consists of gravel and dark-yellow grayish sand from fine to coarse-grained, with rounded and subrounded clasts of lapilli, crystals and weathered pumices. It is a marker of a regressive trend, which permitted coastal progradation and development of river meandering dynamics. Furthermore, in the eastern sector, at a mean depth of 12–14 m b.s.l., a deposit constituted by organic-rich grey mud and peat, which is representative of a backdune or meander cut-off environment, has been identified. The lower degree of permeability of such a deposit characterizes the Marshy Complex (MC). The emplacement of the pyroclastic fall deposits, dating back to the Vesuvius eruption of AD 79, has advanced the aggradation/progradation of the plain. They are composed of both white and gray pumices (Luongo et al. 2003) and are recognizable at depths varying from 8 to DOI 10.1007/s10040-013-1087-8 695 Fig. 3 Key stratigraphic correlation relating lithofacies to hydrofacies for the log B5 in Fig. 2. (a) Borehole stratigraphy; (b) stratigraphic units (Pescatore et al. 2001); (c) hydrostratigraphic-log 15 m b.s.l. along the river pattern and the pre-existing topography. This Pumices Complex (PC), about 2 m thick, is a well-sorted gravel layer, characterized by the highest permeability in the stratigraphic sequence. After the fall event, a pyroclastic density current caused the deposition of ashes from medium-to-fine grained with accretionary lapilli (volcanic sands) that form massive (flow) or laminated structures (surge) (Luongo et al. 2003). These pyroclastic flows and surge deposits allow one to detect two hydrogeological complexes: the Volcanic Ash Complex (VAC) and the Silty Sandy Complex (SSC). The latter is interbedded between two volcanic ash layers, which are at depths of 12 and 7 m b.s.l., respectively. Post AD 79 eruption deposits, representative of mixed pyroclastic and alluvial materials, close the sedimentary succession. They consist of ashes, sands and gravels and Table 1 Order of magnitude of the hydraulic conductivity, K, for the hydrofacies recognized in the investigated subsurface ID Hydrofacies SCC RC USGC Silty Clayey Complex Reworked Complex Upper Sandy Gravelly Complex Volcanic Ash Complex Silty Sandy Complex Pumices Complex Marshy Complex Sandy Gravelly Complex Marine Sands Complex VAC SSC PC MC SGC MSC K (m/s) Aquitard Aquifer/aquitard Aquifer 10–7 10–3–10–5 10–3 Aquiclude Aquitard Aquifer Aquiclude Aquifer 10–9 10–5 10–1 10–9 10–2–10–3 Aquifer 10–2 Hydrogeology Journal (2014) 22: 691–703 define the following hydrofacies: Silty Clayey Complex (SCC), Reworked Silt-Sand-Gravel Complex (RC) and Upper Sandy Gravelly Complex (SGC). The stratigraphic and geomorphologic setting of Quaternary deposits is transferred in a 3D aquifer model shown in Fig. 5, whose hydrostratigraphy points out the distribution of alluvial, marine and volcanic units and the influence of the variation of the coastal onlap. In particular, the representative volume of such complex aquifer system is shown in Fig. 5a with a cut along the 1– 2 transect, which evidences the paleoriverbed, also visible in the hydrostratigraphic cross-section along the BB′ line (Fig. 5c). The section AA′ (Fig. 5b), instead, shows the lateral facies change toward the marine deposits. Figure 5 suggests that, in the eastern area, the groundwater flow is confined and takes place mainly in the deeper and most permeable aquifer, i.e. the SandyGravelly and the Pumices complexes, underlying the volcanic ash layer at depth of about 14 m b.s.l. The top of this aquifer declines toward the NW due to the path of the buried paleoriverbed (Fig. 5a,c), which is revealed by coarser and most permeable reworked sediments. A shallower aquifer is present at depth, variable from 4 to 7 m b.s.l., and is represented by the Upper Sandy Gravelly and the Reworked Silt-Sand-Gravel complexes. The poor grain-size sorting and the discontinuity of the most permeable layers cause this aquifer to be less productive. In the western sector, the aquifer becomes unconfined and unique, because the aforementioned most permeable complexes are interbedded with the Marine Sands Complex (Fig. 5b). In addition, it is worth noting that the buried path DOI 10.1007/s10040-013-1087-8 696 Fig. 4 Hydrogeological fence diagram and main geological information. The codes of the complexes are explained in Table 1 Fig. 5 a 3D hydrofacies model with a cut along the 1–2 line that shows the paleoriverbed. Key hydrostratigraphic cross-sections along the b AA’ and c BB’ lines, showing, respectively, the lateral facies change toward the marine deposits and the buried paleoriverbed Hydrogeology Journal (2014) 22: 691–703 DOI 10.1007/s10040-013-1087-8 697 Fig. 6 a 2D inversion of the resistivity tomography data acquired along the cyan profile shown in the plan-view map of Fig. 5. b 3D inversion of the ERT data acquired along the profiles of Fig. 5 and shown with two cutting planes, respectively, parallel and perpendicular to the N135W oriented profiles of the paleoriverbed, which significantly affects the hydrostratigraphical sequence, constrains the aquifer/aquitard distribution in the NE sector of the survey area, where log data are missing (see inset of Fig. 5). Thus, to improve the paleoriverbed buried pattern and to investigate the spatial continuity of hydrofacies, which is crucial for management and protection of water resources of this highly urbanized costal area, a 2D electrical resistivity tomography (ERT) survey in the northern part of the test site (see Fig. 5) was carried out, as discussed in the following. 3D geoelectrical model In the last few decades, increasing attention has been devoted to applications of geophysical methods to hydrogeological studies. Hydrogeophysics, in fact, is a new branch of applied geophysics (Hubbard and Rubin 2002; Rubin and Hubbard 2005; Vereecken et al. 2006; Hinnel et al. 2010) that is founded on non-invasive geophysical measurements for deriving physical properties which can be correlated to hydrological variables (e.g., water content, porosity, solute concentration). Indeed, conventional sampling techniques, besides being invasive, are generally applied at few measurement points, thus not permitting an exhaustive description of the shallow and deep vadose zones. Conversely, the geophysical techniques provide high-resolution images of hydrogeological Hydrogeology Journal (2014) 22: 691–703 structures in terms of geometrical and physical quantities also in the vadose zone (Daily et al. 1992; Cassiani and Binley 2005). Among all suitable methods, electrical geophysical prospecting is doubtless one of the most prominent and used for hydrogeological purposes, by virtue of the well-known dependence of the electrical resistivity on hydraulic properties (Slater 2007; Di Maio and Piegari 2011). In addition, the development of increasingly effective 2D or 3D data acquisition techniques has allowed detailed space mapping of resistivity values that, if monitored over time, can well describe the spatial and temporal variability of an aquifer system (Sandberg et al. 2002; Vanderborght et al. 2005; Ogilvy et al. 2009). In particular, the geoelectrical prospecting used in the present study was aimed at constraining the hydrofacies model in the upper part of the NE sector of the survey area (see Fig. 5), in order to recognize also in that zone the heterogeneous soil profile with multi-layer structures (see previous section), and to expand the in-depth modelling. Accordingly, a 2D resistivity tomography survey was carried out in an area of about 200×120 m2, where geological investigations were not allowed. Specifically, 2D electrical resistivity tomographies were performed along 12 profiles, five with orientation N45W and 15 m distant from each other, and seven with orientation N135W and 30 m distant from each other (Fig. 5). Profiles that are N45W oriented have a length varying between 68 DOI 10.1007/s10040-013-1087-8 698 Fig. 7 a 2D inversion of the resistivity tomography data acquired along the orange profile shown in the plan-view map of Fig. 5. b 3D inversion of the ERT data acquired along the profiles of Fig. 5 and shown with three different cuts and 115 m; profiles that are N135W oriented have a length varying between 140 and 185 m. The measurements were realized by using an axial dipole-dipole electrode configuration with an inter-electrode distance of 2 or 5 m, in order to attain a good horizontal resolution and data coverage up to the maximum exploration depth (Loke 2004). The apparent resistivity data inversion was performed by using the Res2D algorithm (Loke and Barker 1996; Loke and Dahlin 2002), which consists of a least-squares deconvolution method based on a linear process between acquired apparent resistivity values as a function of real resistivities. Then, the 2D ERT data were also inverted by using the Res3Dinv algorithm (Loke 2004) to obtain a 3D resistivity-data distribution of the investigated volume. Figure 6 shows the inversion results of the 2D ERT measurements carried out along a N135W oriented profile (Fig. 6a) and two clipped volumes of the 3D data- inversion results (Fig. 6b). The upper cross section of Fig. 6b, which corresponds to the selected profile, well matches the 2D image of Fig. 6a. From the figures, a conductive level sandwiched between two relatively highresistivity layers appears in the depth range of about 2.5– 4 m b.g.l. By considering the values that characterize this conductive layer, it can be correlated to the less productive shallow aquifer system identified by the hydrogeological analysis, which is composed of the Silty Clayey and the Upper Sandy Gravelly Complexes, and underlies a quite dry cover soil. In fact, the non-continuity of the conductive layer (see Fig. 6b) well describes the heterogeneity emerging by the borehole log stratigraphy; furthermore, a decrease of the resistivity values clearly appears at the maximum investigation depth. Figure 7 shows the inversion results of the 2D ERT measurements carried out along a N45W oriented profile Fig. 8 Mean resistivity values and variation ranges for a three N45W-oriented profiles and b five N135W-oriented profiles. Black arrows indicate associations with hydrofacies (see also Table 2 and Fig. 9) Hydrogeology Journal (2014) 22: 691–703 DOI 10.1007/s10040-013-1087-8 699 Table 2 Order of magnitude of the hydraulic conductivity and resistivity ranges of the six recognized hydrostratigraphic units Hydrofacies Depth (m) Hydraulic conductivity (m/s) Resistivity (ohm m) Silty Clayey Complex Upper Sandy Gravelly Complex Volcanic Ash Complex, Silty Sandy Complex Pumices Complex Sandy Gravelly Complex 1–5 5–10 10–7 10–3 20–75 75–20 10–15 10–9–10–5 40–10 15 >15 10–1 10–2–10–3 5–20 20–80 (Fig. 7a) and three different cuts for the whole inverted data volume (Fig. 7b). The images show a relatively resistive layer approximately in the first 10 m b.g.l., and a conductive layer with resistivity values of the order of 10 ohm m, which can be ascribed to the presence of an aquifer. Moreover, an increase of the resistivity values approximately starting from 25 m b.g.l. could indicate the transition to a coarser grain size. It is worth noticing that the inter-electrode distance used for N45W oriented profiles, i.e. 5 m, does not allow a detailed definition of the geological strata sequence in the first 10 m of the section. An integrated interpretation of hydrogeological and geophysical data The combination of the geophysical measurements and the hydrogeological analysis showed in previous sections allows for the characterization of the electrical behavior of the recognized hydrofacies and for extension of the hydrofacies model toward the NE sector of the survey area, where hydrogeological information is lacking. As concerns the electrical characterization of the hydrofacies, the variation range of the resistivity values observed for each investigation depth was calculated. As example, Fig. 8 shows the mean values and variation ranges of electrical resistivity for three of the N45W oriented profiles (Fig. 8a) and five of the N135W-oriented profiles (Fig. 8b). Looking at the behaviour of the resistivity mean values with the depth shown in Fig. 8a, it clearly emerges that there is a minimum ranging from 5 to 10 ohm m, identified at about 15 m b.s.l. and representative of the electrofacies of the Pumice Complex. In addition, the increasing resistivity trend at deeper depths permits one to locally assess the permeability changes due to the grain-size sorting within the Sandy Gravelly Complex. Namely, it maps the coarse sand and gravel layer, whose thickness is at least 10 m. Such a feature allows for expansion of the hydrofacies model in both lateral and vertical directions. Conversely, the shallower portion of the hydrostratigraphic succession, about 2.5 m thick and composed of the Silty Clayey and the Upper Sandy Gravelly complexes, shows a relatively strong variation of the resistivity values ranging from about 20 to 100 ohm m. This variability is likely attributable to a different saturation degree of the aquifer system due to both the hydrofacies grain-size sorting and the aquifer and aquitard intercalation. As a result, it is expected that the Silty Clayey and the Upper Sandy Fig. 9 Correlation between geoelectrical and hydrogeological layers. a hydrostratigraphic cross-section along the BB’ line of Fig. 5. b 3D resistivity model with a cut plane parallel to the N135W oriented profiles Hydrogeology Journal (2014) 22: 691–703 DOI 10.1007/s10040-013-1087-8 700 Fig. 10 3D hydrogeological–geophysical model along the 1–2 cut of Fig. 5 Gravelly complexes exhibit permeability values less than that of the deeper layers. In the first 2 m of depth, Fig. 8b shows a resistivity decreasing trend (from about 400 to 500 ohm m) corresponding to Infill and Soil, and a slight increase of the resistivity mean values down to about 5 m b.s.l., ascribable to the Silty Clayey Complex. Interestingly, such an increasing trend is also visible in Fig. 8a, even if it is stretched due to the reduced resolution related to the resistivity data-sampling step. Finally, from about 5 to 10 m b.s.l., the resistivity decreases due to the presence of the Upper Sandy Gravelly Complex. The knowledge of porosity and water salinity values of the recognized hydrogeological complexes would strengthen the proposed correlation between electrolayers and hydrofacies. However, such values are not available and invasive tests cannot be performed in the survey area. Besides, laboratory and in situ measurements of these parameters are available for the Sandy Gravelly Complex in neighbouring sites. In particular, a water electrical conductivity value of about 1,400 μS/cm and an average porosity of 53 % have been retrieved. By using such a conductivity value and the average resistivity value of 50 Fig. 11 Comparison between hydrostratigraphic setting and groundwater conditions Hydrogeology Journal (2014) 22: 691–703 DOI 10.1007/s10040-013-1087-8 701 ohm m (see Table 2), an average porosity value of 58 % is estimated through the empirical Archie’s law (Archie 1942). This value is in very good agreement with the observed experimental value. To summarize the results discussed previously, Fig. 9 relates the results coming from hydrogeological and geophysical studies and Table 2 provides a summary of the outlined features. Finally, for the stratigraphic framework derived from the 2D electrical resistivity survey in the NE sector, an additional inversion of hydrogeological data was performed. As a result, the hydrogeological model was improved and the pathway of the paleoriverbed was clearly defined (Fig. 10). To validate such a hydrogeological–geophysical model and to trial the hydrodynamic conditions of the survey area, groundwater level measurements in 20 piezometers were carried out. They have shown, in the eastern sector, the shallower aquifer is dry because of widespread overpumping, and groundwater flow can take place only after significant rainy events. Furthermore, groundwater flow is detected within the deeper Sandy Gravelly Complex in the eastern sector, and within the Marine Sands Complex in the western sector. The water table map (Fig. 11) shows that groundwater flows into the sea along a main NE–SW direction. Conversely, it is noteworthy to observe the influence of aquifer system anisotropy on groundwater flowpaths. According to the 3D hydrofacies model, it is possible to correlate the main flowpath to coarser sediments and most permeable hydrofacies along the buried paleoriver channel, as well as to associate the groundwater divide with the recognized Marshy Complex aquiclude at mean depth of 13 m b.s.l. (Fig. 11). Conclusions In this report, the potential of an integrated and interdisciplinary approach to the assessment of complex multilayered aquifer systems is explored. Geophysical and hydrogeological data coming from close survey areas have been combined to identify the architecture of a complex system made of five electrolayers, each of them corresponding to one or more hydrofacies. In particular, the Silty Clayey and the Upper Sandy Gravelly complexes, which are expected to exhibit permeability values less than that of the deeper layers; the Pumice Complex, characterized by the lowest values of resistivity ranging from 5 to 10 ohm m at a depth of about 15 m; a Sandy Gravelly Complex, whose increasing trend of resistivity values permits local assessment of the permeability changes due to the grain-size sorting with the depth. As concerns the groundwater flow system, a main flow was detected within the deeper Sandy Gravelly Complex in the eastern sector, and within the Marine Sands Complex in the western sector, flowing into the sea along a main NE– SW direction. It is worth noting that the realistic reconstruction of a 3D hydrogeophysical model has provided continuous Hydrogeology Journal (2014) 22: 691–703 information in space, beyond the availability of logs, in a highly urbanized area, where measurements are very often prevented by logistic conditions, and management and protection of water resources are strongly required. Acknowledgements Authors are grateful to Pasquale Paduano for his very friendly collaboration. Authors wish to thank Roberto Balia and an anonymous reviewer for their useful suggestions and comments, which helped improve the manuscript. References Ahmed S, Marsily G, Talbot A (1988) Combined use of hydraulic and electrical properties of an aquifer in a geostatistical estimation of transmissivity. Groundwater 26:78–86 Aprile F, Toccaceli R, Sbrana A (2004) Il ruolo dei depositi piroclastici nell’analisi cronostratigrafica dei terreni quaternari del sottosuolo della Piana Campana [The role of pyroclastic deposits in the chronostratigraphic analysis of Quaternary terrains in the Campania Plain subsoil]. Il Quaternario 17:547– 554 Archie GE (1942) The electrical resistivity log as an aid in determining some reservoir characteristics. Trans AIME 146:54–62 Balducci S, Vaselli M, Verdini G (1983) Exploration well in the Ottaviano permit, Italy, Trecase 1. Eur. Geot. Update, 3rd Int. Sem., Munich, 29 November–1 December 1983, pp 407 – 418 Barberi F, Innocenti F, Lirer L, Munno R, Pescatore T, Santacroce R (1978) The Campanian Ignimbrite: a major prehistoric eruption in the Neapolitan area (Italy). Bull Volcanol 41:1–22 Barberi F, Cioppi D, Ghelardoni R, Nannini R, Sommaruga C, Verdiani G (1980) Integrated geothermal reconaissance of the Somma-Vesuvius system. 2nd Int. Sem. Results E.C. Geothermal Energy Res., Strasbourg, March 1980, pp 141–153 Barra D, Romano P, Santo A, Campaiola L, Roca V, Tuniz C (1996) The Versilian transgression in the Volturno river plain (Campania, southern Italy): paleoenvironmental history and chronological data. Il Quaternario 9:445–458 Beven K (2007) Towards integrated environmental models of everywhere: uncertainty, data and modeling as a learning process. Hydrol Earth Syst Sci 11:460–467 Brancaccio L, Cinque A, Romano P, Rosskopf C, Russo F, Santangelo N, Santo A (1991) Geomorphology and neotectonic evolution of a sector of the Tyrrhenian flank of the southern Apennines (region of Naples, Italy). Zeitschrift Geomorphol Suppl 82:47–58 Brocchini D, Principe C, Castratori D, Laurenzi MA, Gorla L (2001) Quaternary evolution of the southern sector of the Campania Plain and early Somma-Vesuvius activity: insights from the Trecase 1 well. Mineral Petrol 73:67–91 Cassiani G, Binley AM (2005) Modeling unsaturated flow in a layered formation under quasi-steady state conditions using geophysical data constraints. Adv Water Resour 28:467–477 Cassiani G, Böhm G, Vesnaver A, Nicolich R (1998) A geostatistical framework for incorporating seismic tomography auxiliary data into hydraulic conductivity estimation. J Hydrol 206:58–74 Catuneanu O (2006) Principles of sequence stratigraphy. Elsevier, Amsterdam Celico F, Piscopo V (1995) Idrodinamica sotterranea e vulnerabilità all’inquinamento delle piane del Sarno e della Solofrana (Campania) [Groundwater flow and pollution vulnerability of the Sarno and Solofrana River plains (Campania Region)]. Quad Geol Appl 2:407–414 Chandra S, Ahmed S, Ram A, Dewandel B (2008) Estimation of hard rock aquifers hydraulic conductivity from geoelectrical measurements: a theoretical development with field application. J Hydrol 357:218–227 DOI 10.1007/s10040-013-1087-8 702 Chen J, Hubbard S, Peterson J, Williams K, Fienen M, Jardine P, Watson D (2006) Development of a joint hydrogeophysical inversion approach and application to a contaminated fractured aquifer. Water Resour Res 42, W06425 Cinque A, Patacca E, Scandone P, Tozzi M (1993) Quaternary kinematic evolution of the Southern Apennines: relationships between surface geological features and deep lithospheric structures. Annal Geofis 36:249–259 Daily WA, Ramirez A, Labrecque D, Nitao J (1992) Electrical resistivity tomography of vadose water movement. Water Resouces Res 28:1429–1442 De Marsily G, Delay F, Goncalvès J, Renard P, Teles V, Violette S (2005) Dealing with spatial heterogeneity. Hydrogeol J 13:161–183 De Vivo B, Rolandi G, Gans PB, Calvert A, Bohrson WA, Spera FJ, Belkin HE (2001) New constraints on the pyroclastic eruptive history of Campanian volcanic plain (Italy). Miner Petrol 73:47–65 Di Maio R, Piegari E (2011) Water storage mapping of pyroclastic covers through resistivity measurements. J Appl Geophys 75:196–202 Di Maio R, Mauriello P, Patella D, Petrillo Z, Piscitelli S, Siniscalchi A (1998) Electric and electromagnetic outline of the Mount Somma-Vesuvius structural setting. J Volcanol Geotherm Res 82:219–238 Drahor MG, Berge MA, Göktürkler G, Kurtulmuş TO (2011) Mapping aquifer geometry using electrical resistivity tomography: a case study from Şanlıurfa, south-eastern Turkey. Near Surf Geophys 9:55–66 Fabbrocino S, Lanari R, Celico P, Termolini G, Zeni G (2007) Groundwater pumping and land subsidence in the Sarno River plain. Mem Descr Carta Geol Ital 76:163–174 Felletti F, Bersezio R, Giudici M (2006) Geostatistical simulation and numerical upscaling to model groundwater flow in a sandygravel, braided river, aquifer analogue. J Sedimen Res 76:1215– 1229 Ferré T, Bentley L, Binley A, Linde N, Kemna A, Singha K, Holliger K, Huisman JA, Minsley B (2009) Critical steps for the continuing advancement of hydrogeophysics. EOS Trans AGU 90:200. doi:10.1029/2009EO230004 Galloway DL (2010) The complex future of hydrogeology. Hydrogeol J 18:807–810 Gasparini P (1998) Looking inside Mt Vesuvius. EOS Trans AGU 79:229–232 Hinnel AC, Ferrè TPA, Vrugt JA, Huisman JA, Moysey S, Rings J, Kowalsky MB (2010) Improved extraction of hydrological information from geophysical data through coupled hydrogeophysical inversion. Water Resour Res. doi:10.1029/ 2008WR007060 Hsien-Tsung L, Yih-Chi T, Chu-Hui C, Hwa-Lung Y, Shih-Ching W, Kai-Yuan K (2010) Estimation of effective hydrogeological parameters in heterogeneous and anisotropic aquifers. J Hydrol 389:57–68 Hubbard S, Rubin Y (2002) Hydrogeophysics: state-of-the-discipline. EOS Trans AGU 83:602–606 Hubbard S, Rubin Y, Majer E (1999) Spatial correlation structure estimation using geophysical and hydrogeological data. Water Resour Res. doi:10.1029/1999WR900040 Huntley D (1986) Relations between permeability and electrical resistivity in granular aquifers. Groundwater 24:466–474 Hyndman DW, Harris JM, Gorelick SM (1994) Coupled seismic and tracer test inversion for aquifer property characterization. Water Resour Res 30:1965–1977 Kelly WE, Frohlich KR (1985) Relation between aquifer electrical and hydraulic properties. Groundwater 23:182–189 Koch K, Wenninger J, Uhlenbrook S, Bonnell M (2009) Joint interpretation of hydrologic and geophysical data: electrical resistivity tomography results from a process hydrological research site in the Black Forest Mountains Germany. Hydrol Process 23:1501–1513 Kowalsky MB, Finsterle S, Peterson J, Hubbard S, Rubin Y, Majer E, Ward A, Gee G (2005) Estimation of field‐scale soil Hydrogeology Journal (2014) 22: 691–703 hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour Res. doi:10.1029/ 2005WR004237 Laney RL, Davidson CB (1986) Aquifer nomenclature guidelines. US Geol Surv Open-File Rep 86-594 La Torre P, Nannini R, Sbrana A (1983) Geothermal exploration in southern Italy: geophysical interpretation of Vesuvian area. Boll Geofis Teor Appl 26:197–208 Lemon AM, Jones NL (2003) Building solid models from boreholes and user defined cross sections. Comput Geosci 29:547–555 Linde N, Binley A, Tryggvason A, Pedersen LB, Revil A (2006) Improved hydrogeophysical characterization using joint inversion of crosshole electrical resistance and ground penetrating radar traveltime data. Water Resour Res. doi:10.1029/2006WR005131 Loke MH (2004) Tutorial: 2-D and 3-D electrical imaging surveys. http://www.cas.umt.edu/geosciences//faculty/sheriff/495-subsurface imaging in archaeology/Sources/Loke_elect_tutorial.pdf. Accessed 25 September 2013 Loke MH, Barker RD (1996) Rapid least-squares inversion of apparent resistivity pseudosections using a quasi-Newton method. Geophys Prospect 44:131–152 Loke MH, Dahlin T (2002) A comparison of the Gauss-Newton and the quasi-Newton methods in resistivity imaging inversion. J Appl Geophys 49:149–162 Luongo G, Perrotta A, Scarpati C (2003) Impact of AD 79 explosive eruption on Pompeii. I. Relationship among depositional mechanisms of the pyroclastic products, the framework of the buildings and the associated destructive events. J Volcanol Geoth Res 126:201–223 Marturano A, Aiello G, Barra D (2011) Evidence for late Pleistocene uplift at the Somma-Vesuvius apron near Pompeii. J Volcanol Geoth Res 202:211–227 Mele M, Bersezio R, Giudici M (2012) Hydrogeophysical imaging of alluvial aquifers: electrostratigraphic units in the quaternary Po alluvial plain (Italy). Int J Earth Sci 101:2005– 2025 Ogilvy RD, Kuras O, Meldrum PI, Wilkinson PB, Chambers JE, Sen M, Ginsbert J, Jorreto S, Frances I, Pulido-Bosh A, Tsourlos P (2009) Automated time-lapse electrical resistivity tomography (ALERT) for monitoring coastal aquifers. Near Surf Geophys 7:367–375 Ortolani F, Aprile F (1985) Principali caratteristiche stratigrafiche e strutturali dei depositi superficiali della Piana Campana [Main stratigraphic and structural features of shallow deposits of the Campania Plain]. Boll Soc Geol Ital 104:195–206 Ouellon T, Lefebvre R, Marcotte D, Boutin A, Blais V, Parent M (2008) Hydraulic conductivity heterogeneity of a local deltaic aquifer system from the kriged 3D distribution of hydrofacies from borehole logs, Valcartier, Canada. J Hydrol 351:71–86 Patacca E, Sartori R, Scandone P (1990) Tyrrenian basin and apenninic arcs: kinematic relations since late Tortonian times. Mem Soc Geol Ital 45:425–451 Pescatore T, Senatore MR, Capretto G, Lerro G (2001) Holocene costal environments near Pompeii before A.D. 79 eruption of mount Vesuvius (Italy). Quat Res 55:77–85 Poland JF, Lofgren EE, Riley FS (1972) Glossry of selected terms useful in studies of the mechanics of aquifer systems and land subsidence due to fluid withdrawal. US Geol Surv Water Suppl Pap 2025, 9 pp Robins NS, Rutter HK, Dumpleton S, Peach DW (2005) The role of 3D visualization as an analytical tool preparatory to numerical modeling. J Hydrol 301:287–295 Rolandi G, Bellucci F, Heizler MT, Belkin HE, De Vivo B (2003) Tectonic controls on the genesis of ignimbrites from the Campanian Volcanic Zone, southern Italy. Miner Petrol 79:3–31 Romano P, Santo A, Voltaggio M (1994) L’evoluzione geomorfologica della pianura del F. Volturno (Campania) durante il tardo Quaternario (Pleistocene medio-superiore Olocene) [Geomorphological evolution of the Volturno River plain (Campania) in the late Quaternary (mid-late Pleistocene Holocene)]. Il Quaternario 7:41–56 DOI 10.1007/s10040-013-1087-8 703 Rubin Y, Hubbard S (2005) Inverse and forward modeling using multiple data sets: the hydrogeophysical challenge. In: Rubin Y, Hubbard S (eds) Hydrogeophysics. Springer, Dordrecht, The Netherlands Sandberg SK, Slater LD, Versteeg R (2002) An integrated geophysical investigation of the hydrogeology of an anisotropic unconfined aquifer. J Hydrol 267:227–243 Santacroce R, Cioni R, Marianelli P, Sbrana A, Sulpizio R, Zanchetta G, Donahue DJ, Joron JL (2008) Age and whole rock-glass compositions of proximal pyroclastics from the major explosive eruptions of Somma-Vesuvius: a review as a tool for distal tephrostratigraphy. J Volcanol Geoth Res 177:1–18 Sinha R, Israil M, Singhal DC (2009) A hydrogeophysical model of the relationship between geoelectric and hydraulic parameters of anisotropic aquifers. Hydrogeol J 17:495–503 Slater L (2007) Near surface electrical characterization of hydraulic conductivity: from petrophysical properties to aquifer geometries—a review. Surv Geophys 28:169–197 Tronicke J, Holliger K (2005) Quantitative integration of hydrogeophysical data: conditional geostatistical simulation for characterizing heterogeneous alluvial aquifers. Geophysics 70:H1–H10 Hydrogeology Journal (2014) 22: 691–703 View publication stats Vanderborght J, Kemna A, Hardelauf H, Vereecken H (2005) Potential of electrical resistivity tomography to infer aquifer transport characteristics from tracer studies: a synthetic case study. Water Res Resour. doi:10.1029/2004WR003774 Vereecken H, Binley A, Cassiani G, Revil A, Titov K (2006) Applied hydrogeophysics. Springer, Dordrecht, The Netherlands. doi:10.1007/978-1-4020-4912-5_1 Vienken T, Dietrich P (2011) Field evaluation of methods for determining hydraulic conductivity from grain size data. J Hydrol 400:58–71 Zollo A, D’Auria L, De Matteis R, Herrero A, Virieux J, Gasparini P (2001) Bayesian estimation of 2-D P-velocity models from active seismic arrival time-data: imaging of the shallow structure of Mt. Vesuvius (southern Italy). Geophys J Int 151:506–582 Zollo A, Gasparini P, Virieux J, Le Meur H, De Natale G, Biella G, Boschi E, Capuano P, De Franco R, Dell’Aversana P, De Matteis R, Guerra I, Iannaccone G, Mirabile L, Vilardo G (1996) Seismic evidence for a low-velocity zone in the upper crust beneath Mount Vesuvius. Science 274:592– 594 DOI 10.1007/s10040-013-1087-8