Subsurface geological and geophysical data from the Po Plain and the northern Adriatic Sea (north Italy)

The Po Plain (Italy) is one of the most densely populated and productive regions of Europe, characterized by a flourishing economy (also linked to strategic subsurface resources) and several world cultural and natural heritage sites. The coupling of socio-economic interests with geological hazards (i.e. seismic, subsidence, and flooding hazards) in this area requires accurate knowledge of the subsurface geology, the active geological processes, and the impact of human activities on natural environments to mitigate the potential natural and anthropic risks.

Most data unveiling the subsurface geology of this region were produced by the hydrocarbon exploration industry. Indeed, the Po Plain hosts many hydrocarbon fields that have been discovered since the early 1950s, giving rise to the subsurface exploration through extensive seismic reflection surveys and drilling of numerous deep wells. In this work, geological and geophysical data from 160 deep wells drilled for hydrocarbon exploration and/or exploitation purposes in the Po Plain and in the facing northern Adriatic Sea have been collected and digitized along with several published geological cross-sections and maps. These data have been used to reconstruct the overall subsurface 3D architecture and to extract the physical properties of the subsurface geological units.

The digitized data are suitable to be imported into geo-software environments so as to derive the geophysical and mechanical properties of the geological units for a wealth of applied and scientific studies such as geomechanical, geophysical, and seismological studies.

The integrated dataset may represent a useful tool in defining regional first-order strategies to ensure the safety of the urbanized areas and human activities and to reduce natural and anthropic risks that may affect this crucial region of Europe. In particular, the data collected would be useful to highlight sensible areas where data collection and more detailed studies are needed. Nowadays, such issues are particularly relevant for the underground industry development related to the increasing interest in possible CO₂ and hydrogen underground storage, which can play a fundamental role in the energy transition process towards decarbonization goals. The full dataset is available at the following link: https://doi.org/10.5281/zenodo.8126519 (Livani et al., 2023).

1 Introduction

The Po Plain in the north of Italy is the most productive and prosperous region of Italy (Fig. 1a), with a per-capita income similar to that of central and northern European countries (https://ec.europa.eu/eurostat, last access: 8 July 2022; Helliwell and Putnam, 1995; Tabellini, 2010; European Commission, 2016). The area is densely populated (https://ec.europa.eu/eurostat, last access: 8 July 2022; Fig. 1a) and hosts numerous UNESCO world heritage sites such as the cities of Venice, Verona, Bergamo, and Ravenna (https://whc.unesco.org/, last access: 13 December 2022; Fig. 1b).

Since the second half of the 20th century, the discovery and exploitation of numerous hydrocarbon (mostly gas) resources contributed to the economic development of Italy. Indeed, the Po Plain and the facing northern Adriatic Sea host the majority of hydrocarbon fields in the country (extracting almost 33 % of the total national gas production) and most of the Italian underground gas storage sites, which have been operative since the 1960s (https://unmig.mise.gov.it/, last access: 9 January 2023; Fig. 1b).

The Eni-AGIP Company hydrocarbon exploration and exploitation activities in the Po Plain and northern Adriatic Sea led to the production of a large amount of subsurface data, which include: (i) seismic data acquired through extensive regional 2D and 3D seismic surveys for the development of onshore and offshore hydrocarbon fields and (ii) well data acquired during the drilling of explorative and developmental wells. The subsurface dataset provided structural, stratigraphic, and sedimentological information, which allowed the obtainment of accurate knowledge of the regional subsurface architecture and geological evolution (e.g. Pieri and Groppi, 1981; Cassano et al., 1986; Casero, 2004; Ghielmi et al., 2010, 2013; Fantoni and Franciosi, 2009, 2010).

The geological evolution of the Po Plain is accompanied by seismic, subsidence, and flooding events that are the main phenomena to consider in natural-hazard assessments. The cultural, social, and economic relevance of this region demands a careful evaluation of both natural hazards and human activities to mitigate the potential correlated risks and to guarantee the safety of the urbanized areas and human activities themselves.

The Po Plain is characterized by moderate seismicity (International Commission on Hydrocarbon Exploration and Seismicity in the Emilia region, 2014). Nonetheless, instrumental and historical intermediate–strong (M_w≥5.0) earthquakes hit the area (Rovida et al., 2022), such as the events that occurred in 1117 (M_w 6.5; Verona area); in 1570 (M_w 5.6; Ferrara area); and 2012, when a seismic sequence in the Modena and Ferrara provinces occurred (main shock M_w 6.09). The latter event promoted the implementation of the seismic monitoring network as an essential instrument for the safe management of industrial activities. Notice that no case of anthropogenically induced seismicity has been documented so far in the study area (Braun et al., 2020).

The Po Plain is also characterized by significant subsidence of both natural and anthropic origin; this became particularly intense after the economic growth following World War II. The natural component is the result of several geological processes, such as sediment load and compaction, vertical tectonic movements, and rebound effects after the last deglaciation (Carminati and Di Donato, 1999). Anthropogenic subsidence is primarily influenced by groundwater withdrawal for industrial, agricultural, and civil uses (e.g. Herrera-Garcia et al., 2021). Subordinately, it is linked to hydrocarbon (mainly gas) exploitation from onshore and offshore reservoirs (Carminati and Martinelli, 2002; Bitelli et al., 2020).

Subsidence mainly affects the central and eastern sectors of the Po Plain: the highest total subsidence rates, greater than −60 mm yr⁻¹, were evaluated for the Bologna city area (Carminati and Martinelli, 2002; Zerbini et al., 2007; Baldi et al., 2009), whereas the highest natural component reaches −2.0 mm yr⁻¹ at the Pede–Apennine zone, near the city of Bologna. The difference between the total (natural plus anthropogenic) vs. natural components of the present-day subsidence shows that the main factor controlling modern subsidence in this region is anthropogenic (Carminati and Martinelli, 2002).

Subsidence is monitored through the implementation of different technologies (GPS, InSAR (interferometric synthetic aperture radar), levelling surveys) that are useful for the prevention, mitigation, and control of the natural processes, as well as of the human ones (e.g. Dacome et al., 2015; Benetatos et al., 2020).

Flood events related to the Po River (which crosses the whole Po Plain) and its tributaries represent an additional natural risk (Domeneghetti et al., 2015, and references therein). In the past century, the worldwide increase of flood susceptibility and risk increased together with the rise of subsidence due to groundwater depletion (e.g. Herrera-Garcia et al., 2021). In the Po Plain, subsidence and flood hazards are presumably linked since there is a clear-cut correlation between high flood frequency and rapid subsidence, whereas only a few floods occurred in low-subsidence areas (Carminati and Martinelli, 2002). Nevertheless, in the first half of 2022, the most significant drought of at least 70 years was observed for the area, causing extensive damage to agriculture and encouraging the entry of the saline wedge at the mouth of the Po River along the Adriatic coast.

The protection of the cultural, social, and economic heritage of the Po Plain makes the adoption of all available measures necessary to reduce the impact of natural and anthropic-derived hazards in the study area.

To this end, the definition of the physical and geometrical attributes of the outcropping and subsurface geological units provides fundamental knowledge for several scientific issues and could be used in a preliminary definition of engineering operations. As an example, the integration of physical parameters of subsurface geological units into a well-defined 3D model can be applied to (i) reduce the uncertainties for earthquake location, contributing to the calculation of more accurate focal mechanisms and performing wave propagation and ground motion simulations (e.g. Magistrale et al., 1996; Süss et al., 2001; Molinari et al., 2015; Livani et al., 2022); and (ii) understand, simulate, and predict the response of the geological body to subsurface natural and anthropic processes. The latter is the case for 3D geomechanical numerical models, which represent effective tools for the evaluation and prediction of the possible effects – both at the surface and in the subsurface – of geofluid extraction and storage to guarantee the safe management of such activities, as well as to quantify and better understand the ongoing geological processes (e.g. tectonic deformation and natural subsidence; Teatini et al., 2006; Codegone et al., 2016; Benetatos et al., 2020). As an example, numerical models can play a fundamental role during the geological sequestration of CO₂ and gas storage (e.g. hydrogen as an energy carrier) in natural underground formations since they are mandatory for the optimization of the development strategies, the maximization of storage efficiency, and monitoring activities.

In this paper, we present an integrated database of geological and geophysical data regarding the subsurface of the Po Plain and the facing Adriatic Sea. The database provides a collection of data distinguished as primitive and derived. The primitive data consist of a detailed and accurate collection and digitization of subsurface information extracted from wells, geological cross-sections, and geological maps. The derived data are obtained from the revision and processing of the primitive data and by gridded surfaces representing the main geological units of the Po Plain subsurface.

It is worth mentioning that, in the study area, the geological literature offering interpretations of the subsoil is abundant in terms of geological reconstructions. Our database represents a collection of the main published works regarding the Po Plain. Detailed studies related to specific sectors of the Po Plain might not be present in our database. Recently, other subsurface geological models have been elaborated upon, providing several isobath maps of the Po Plain (e.g. Turrini et al., 2014; ISPRA, 2015; Molinari et al., 2015; Amadori et al., 2019; D'Ambrogi et al., 2023). It is worth mentioning two projects: the recent GO-PEG project (D'Ambrogi et al., 2023), which provides the geometry of four stratigraphic horizons in the Po Plain area deriving from the GeoMol project (http://www.geomol.eu, last access: 30 June 2023; ISPRA, 2015; Maesano and D'Ambrogi, 2016) and the GeoERA-HotLime project (https://geoera.eu/projects/hotlime6/, last access: 30 June 2023), and the Mambo project (Molinari et al., 2015), which covers almost the entire Po Plain. In our work, due to technical difficulties (e.g. the impossibility of precisely defining the depth reported in the maps), we were not able to integrate all the previous datasets. Anyhow, where possible, we integrated the available regional geological models with other primitive public data to coherently define the geometry of five surfaces of the Po Plain subsurface, which are, from the oldest to the youngest, as follows: the magnetic basement top, the top of the carbonate succession, the Pliocene base, the Calabrian base, and the base of recent continental deposits. The above-mentioned surfaces represent lithological boundaries rather than chronostratigraphic or formational limits, defining units that are expected to show different mechanical behaviour. Most importantly, the gridded surfaces of our dataset come with a series of detailed geophysical and geological parameters extracted from the composite logs of deep wells (i.e. the primitive data). In addition to the geophysical and geological data, a statistical analysis is reported and discussed as a preliminary elaboration of the collected data.

We believe that our dataset provides an important contribution to a broad audience of policymakers and scientists for understanding and evaluating geological and anthropic processes in the area, for setting secure development strategies for correct territory management, and for the reduction of social and economic risks in a strategic area for Italy and Europe.

Figure 2 Study area and structural geology of the Po Plain. (a) Simplified structural map of the Po Plain (modified after Livani et al., 2018). The main buried thrusts of the northern Apennines (red colour) and southern Alps (blue colour) fronts are shown along with the instrumental and historical seismicity (ISIDe Working Group, 2010; Rovida et al., 2022). The white polygon represents the study area. (b) Schematic geological cross-section across the Po Plain along a N–S-oriented line (by Livani et al., 2018). The northern Apennines (on the left) and southern Alps (on the right) fronts can be identified under the Plio–Pleistocene sedimentary cover (in yellow) filling the Po Plain foredeep. Section trace is reported with a dashed red line in panel (a).

2 Geological setting

The Po Plain and the facing Adriatic Sea lie on the buried sector of the Adria microplate, a promontory of the Africa plate or an independent microplate, interposed between the NE-verging northern Apennines and the S-verging southern Alps (Fig. 2; Dercourt et al., 1986; Pieri and Groppi, 1981; Castellarin et al., 1985; Castellarin and Vai, 1986; Doglioni, 1993; Carminati et al., 2003; Carminati and Doglioni, 2012; Fantoni and Franciosi, 2009, 2010; Turrini et al., 2016; Pezzo et al., 2020). The development of these two facing fold-and-thrust belts, connected with the broad collision of the Eurasian and African plates, led to the formation of the Po Plain basin representing the foreland basin of both orogens. Compressional tectonics have affected the area since the middle Eocene time, with the development of WNW–ESE-oriented thrusts in the southern Alps followed, from the Oligocene–lower Miocene onward, by the NW–SE-oriented thrust system of the northern Apennines (Coward and Dietrich, 1989; Carminati and Doglioni, 2012; Carminati et al., 2012; Maesano and D'Ambrogi, 2016).

The structural and sedimentary framework of the area has been constrained using numerous seismic reflection profiles and deep-well logs (e.g. Pieri, 1983; Cassano et al., 1986; Fantoni and Franciosi, 2010; Turrini et al., 2014; ISPRA, 2015; Livani et al., 2018; Amadori et al., 2019; D'Ambrogi et al., 2023, and references therein).

Figure 3 Synthetic stratigraphic columns of the Po Plain. From left to right: the stratigraphy of Emilia and the southern Alps, Parma-Mantova, and Romagna coastal areas is shown (see Fig. 2a). Sequences of the carbonate passive margin (from violet to green colours) and silicoclastic active margin (from brown to light-yellow colours) are distinguished. The main décollement levels are highlighted (red for basal décollement and orange for shallow décollement). The stratigraphic column to the left shows the different fillings and the different formation thicknesses of the southern Alps and the northern Apennines foredeep.

The Apennines front is characterized by three main orogenic arcs from west to east: Monferrato, Emilia, and Ferrara arcs (Fig. 2a). The southern Alps represent the non-metamorphic retrobelt of the double-verging Alpine chain, and, on the western side of the study area, it reaches the southernmost extent with its edge very close to the northern Apennines front (Ravaglia et al., 2006; Fantoni and Franciosi, 2010; Toscani et al., 2014; Fig. 2b). The external fronts of the two facing chains are mostly buried under a siliciclastic sequence (late Eocene–actual) consisting of syntectonic sediments and recent alluvial sediments of the Po River (Pieri and Groppi, 1981; Boccaletti et al., 1985; Bigi et al., 1990a, b; Fantoni and Franciosi, 2010; Ghielmi et al., 2010, 2013; Carminati and Doglioni, 2012; Amadori et al., 2019, and references therein). In detail, the siliciclastic sequence (Fig. 3) can be subdivided into a lower (late Eocene–early Messinian) and an upper (late Messinian–present) cycle (Ricci Lucchi, 1986). The lower cycle, primarily fed by the Alpine chain, consists of silty and shaly deposits (i.e. Gallare Marls; late Eocene to Miocene time) in places intercalated by or interdigitated with sandy and conglomeratic deposits (i.e. Gonfolite Fm.; Oligocene time), passing upward to sandy marls (Marnoso-arenacea Fm.; Langhian to Messinian time), clays (i.e. Colombacci fm.; Messinian time), and evaporitic deposits (i.e. Gessoso Solfifera fm.; Messinian time). The upper cycle, fed by both the northern Apennines and the southern Alps, mainly consists of marine sandy and conglomeratic formations (e.g. Sergnano Gravel, Porto Corsini, Porto Garibaldi, Santerno, and Asti sandstones; Pliocene to middle–late Pleistocene time) and alluvial deposits (middle–late Pleistocene to present; e.g. Muttoni et al., 2003; Garzanti et al., 2011; Ghielmi et al., 2010, 2013; Livani et al., 2018). South of the Po River, the continental deposits consist of alluvial fan and plain deposits embedded in clays and showing elongated shapes, whereas, to the north of the Po River, the sedimentary bodies are wider, generally tabular, and with minor amounts of fine-grained sediments (Ori, 1993; Amorosi and Milli, 2001). These clastic sequences are superimposed on a carbonate and marly substratum (Triassic–middle Eocene), which lies on top of the platform and continental Permian–Triassic formations, lying in turn on the Variscan crystalline basement (for more details about the carbonate sequence, refer to Livani et al., 2018, and references therein). The Triassic deposits are sometimes interposed by intra-sedimentary volcanic bodies (e.g. Pieri and Groppi, 1981; Castellarin, 1985; Cassano et al., 1986; Ghielmi et al., 2010; Livani et al., 2018). The southern Alps buried front consists of a repetition of Cenozoic clastic units stacked above a regional detachment in the Marne di Gallare Fm. (late Eocene) and followed at depth by deep thrust cutting the Mesozoic carbonates and the Variscan crystalline basement (Fantoni et al., 2004; Figs. 2b and 3). The northern Apennine chain developed by off-scraping the Meso–Cenozoic sedimentary cover of the subducting Adria plate made up of siliciclastic, evaporitic, and shallow- to deep-water carbonate deposits (Cati et al., 1987; Bertotti et al., 1993; Casero et al., 1990; Zappaterra, 1990; Grandić et al., 2002; Fantoni and Scotti, 2003; Fantoni and Franciosi, 2010; Masetti et al., 2012; Fig. 3).

Due to the Apennine compressional deformation, the northern Apennines thrust system migrated toward the foreland over time (among many others, see Malinverno and Ryan, 1986; Doglioni, 1991; Patacca et al., 1990; Royden, 1988; Faccenna et al., 2003; Rosenbaum and Lister, 2004; Scrocca et al., 2006, 2007; Carminati and Doglioni, 2012, and references therein). The convergence is still active, as indicated by the moderate-to-high seismicity which historically affects the Po Plain (maximum M_w between 5.5 and 6.5; Rovida et al., 2020, 2022). Studies on the active stress field in Italy (Montone et al., 2004; Devoti et al., 2008; Cuffaro et al., 2010; Montone et al., 2012), based on the analysis of earthquake focal mechanisms, GPS records, and borehole breakout data, identify active compression in the shallow portion of the Po Plain, with the maximum shortening axis being orthogonal to the main orogenic structures (i.e. NNE–SSW). Locally, the young land morphologies, the presence of different orders of fluvial terraces, and the deviation of some rivers (including the Po River) near the buried active tectonic structures are the tangible evidence of the recent tectonic activity (e.g. Burrato et al., 1999, 2003; Boccaletti et al., 2004a, b; Wilson et al., 2009; Livio et al., 2009; Zuffetti and Bersezio, 2020; Bresciani and Perotti, 2014).

Figure 4 Primitive data collected from public sources (i.e. databases and articles from the literature). (a) Geographical locations of the 160 wells from which well data were collected and digitized. Some areas of the Po Plain are characterized by a low well density, hence the successive primitive data collection (i.e. geological cross-sections and maps) was focused within a slightly reduced area indicated by the white polygon. (b) Traces of the geological cross-sections collected from the literature (the number on each section corresponds to a specific geological cross-section in Table 1 where the data source is specified). (c) Primitive data from geological subsurface maps. The source is reported in the legend.

Table 1 List and type of data collected for each lithological boundary.

Table 2 List of the geological cross-sections collected and used to build the 3D geological model. “Section number” refers to section traces reported in Fig. 4. “Figure number” refers to the figure in the original work. “Section” refers to the number of the section in the corresponding figure of the original work. “Repositioned” indicates whether the section trace has been modified according to the intersections with data located more accurately. “Section length” represents the length in kilometres of each geological sections. The total length of the collected geological cross-sections is 6341 km.

Table 3 Wells collected in the database with their location (x coordinate and y coordinate) and the rotary table elevation (m a.s.l.). The coordinates of wellhead locations are reported in the geographical system used in the database (WGS 84/UTM zone 32N; EPSG: 32632). UOI: well identification number. The deepest lithology unit reached by each well corresponds to the following: (1) recent continental deposits, (2) late Pliocene–Pleistocene deposits, (3) late Miocene–late Pliocene deposits, (4) early–late Miocene deposits, (5) Triassic–Eocene carbonate units, and (6) crystalline basement.

3 Database description

We realized our database by collecting, revising, and digitizing geological and geophysical data, originally in raster format, derived from public sources (i.e. databases and literature works). We collected 160 deep-well composite logs, 61 published geological cross-sections, and 10 geological maps (Fig. 4; Tables 1–3). The data were georeferenced to a common geographical system (WGS 84/UTM zone 32N; EPSG: 32632) by using the open-source QGis software (version 3.16.7; http://www.qgis.org, last access: 15 May 2021) and were then uploaded into a 3D geological modelling software (Petrel^® Software, version 2016.2). We organized the database into two groups, namely primitive data and derived data, containing further hierarchical subdivisions created for better organization and comprehension of the database.

The primitive data contain the product of the digitization of the collected isobath maps, the geological cross-sections, and the well data realized and acquired over a long period and for different purposes. The well data include specific sets of well logs aimed at the geological and/or mechanical characterization of the geological units. Therefore, to achieve the integration of the different sets of data, we analysed them with a focus on the stratigraphy and geological age of the interpreted units. Notice that, in 2009, the Executive Committee of the International Union of Geological Sciences (IUGS) ratified a new subdivision of the Quaternary Period and the Pleistocene Epoch, lowering the age of their base from the top of the Gelasian (1.8 Ma) to its base (2.58 Ma; Gibbard et al., 2010). Most of the collected data (Fig. 4; Tables 1, 2) refer to the pre-2009 chronostratigraphic subdivision including the Gelasian age in the upper Pliocene (Rio et al., 1998). On the contrary, some recent works reinterpreted the sedimentary sequence in the Po Plain and the nearby northern Adriatic Sea (e.g. Ghielmi et al., 2010, 2013; ISPRA, 2015; Maesano and D'Ambrogi, 2016; Amadori et al., 2019) using the new chronostratigraphic subdivision. In our database, since the collected primitive data (Fig. 4; Tables 1, 2) are both prior and successive to the 2009 chronostratigraphic subdivision, we homogenized the data considering the pre-2009 Pleistocene base as the Calabrian base.

We performed an accuracy analysis on the primitive data that unravelled several discrepancies in the interpretation of the subsurface geological horizons. Starting from these observations and considering the well data as the best constraint, we filtered the primitive isobath maps and the geological cross-sections to obtain a coherent dataset. The derived data consist of these filtered isobath maps and the geological cross-sections plus a series of regional surfaces of the main geological units of the Po Plain subsurface. These surfaces were generated without considering the fault occurrences and/or displacements, and they were gridded by means of interpolation of filtered primitive data.

All data (both primitive and derived) are provided in a delimited text file format and are organized according to the data type (i.e. well, geological cross-section, map, or gridded surface).

During the processes of data collection and revision, some published subsurface data reconstructions might not have been included in the database for technical reasons, such as isobath maps being published at an inappropriate resolution, depth contour lines not being properly stated, or geological cross-sections having a different lithostratigraphic and structural scheme compared to the one used in this work.

Figure 5 summarizes the technical procedures and the workflow used to process the data and to develop the database. In the following paragraphs, the methods and the produced primitive data are explained in detail.

4 Primitive data: methods and results

Primitive data derive from the digitization of public data that have been graphically and spatially checked to eliminate errors due to low graphical quality and distortions, scale errors, or bad positioning. We digitized the principal horizons reported in the geological cross-sections; the isobaths of the main geological surfaces represented in the geological maps; and the well locations, comprised of their trajectory along depth, lithological and stratigraphical information, and geophysical logs from well composite logs. We performed the entire workflow (data collection, image georeferencing, data quality check, integration, and model building) by using QGis and Petrel^® software. The digitized data were then organized into text files. Further details about the data processing are given below.

4.1 Well data

The source of well data is the ViDEPI project database (http://www.videpi.com, last access: 1 July 2023). We collected data from 160 well logs, originally in a raster format (scale 1:1000), drilled in the Po Plain and in the northern Adriatic Sea (Fig. 4). Borehole information, such as wellhead coordinates, rotary table elevation, measured depths, true depths, total depth, and deviation survey, is indicated on the composite logs, along with lithological, stratigraphic, and fluid saturation information (Fig. 6). In addition, composite logs include the spontaneous potential log (SP), which is used for lithological characterization and stratigraphic correlation purposes, and the resistivity log (Res), which is used for the identification of hydrocarbon-bearing intervals; in the more recent well master logs, the SP log is replaced by, or in certain cases complemented by, a gamma ray log (GR). Furthermore, 133 out of 160 well composite logs of our dataset also include sonic log registrations that provide insights into sonic velocity variations with depth (Fig. 6). Further lithological information derives either from drill cuttings or from laboratory analysis of core samples collected from wells.

An accurate revision of the wellhead positioning and the subsurface trajectory was necessary to properly collect and organize the well data and to furnish a dataset suitable to be imported into the most common 3D modelling software. From each well, we first transformed the reported geographic coordinates (expressed in ROMA 40 as geodetic datum) into projected coordinates using the WEST BOAGA projection (geodetic datum: ROMA 40) and then converted into the geographical system used for our database (i.e. WGS 84/UTM zone 32N; EPSG: 32632). Most of the wells are vertical; hence, the well trajectory at depth is set using the wellhead coordinates and the total depth values. On the contrary, the path of the directional wells was reconstructed by using inclination (Inc) and azimuth (Az) information at the depth reported in the composite logs.

Table 4 List of collected wells and the relative log availability. “UOI” indicates the identification number for each well. GR: gamma ray log. SP: spontaneous potential log.

We then digitized the available SP log, GR log, and sonic log of each well using the WebPlotDigitizer software (Rohatgi, 2014). Table 4 shows the log availability for each well in the project. The digitization procedure was performed manually with a variable sampling step or by means of a semi-automatic method of line recognition. The digitized logs were then resampled to a constant step of 0.5 m. Figure 7 shows an example of the digitization process.

Figure 5 Workflow used for the database creation.

Figure 6 Typical composite log (scale 1:1000) available from the ViDEPI database (http://www.videpi.com, last access: 1 July 2023), originally in .pdf format. (a) Section of the well name, well coordinates information, and well log legend. (b) Section of the well log data, lithological cuttings, and completion information. (c) Notes section (e.g. well trajectory, core information, geological information, technical data).

Figure 7 Example of the digitization process for well Correggio 39 dirA. Panel (a) shows the original spontaneous potential log in the composite log (raster format) and the digitalized points (in red). Panel (b) shows the digitized (1) SP and (2) sonic logs and the classification of (3) lithology and (4) hydrocarbon-bearing sections.

The resistivity log (Res) was used to identify the mineralized intervals of wells with gas-bearing layers and then to create a new discrete type of log with indications regarding hydrocarbon-bearing and water-bearing layers. Subsequently, Res was used to isolate and remove the sonic log measurements for those geological intervals affected by the hydrocarbon presence (Figs. 7 and 8).

Figure 8 (a) Variation of transit time with depth for all samples showing the effect of the presence of gas. (b–e) Variation of transit time with depth for the sand, marl, shale, and limestone lithologies.

The stratigraphic information on the composite logs together with the SP and GR logs was used to perform a correlation at the regional scale by identifying those units showing different lithological properties and defining surfaces dividing geological successions with different mechanical properties. The main units recognized are outlined (from youngest to oldest) as follows:

• recent clastic deposits of the Po Plain and the Adriatic Sea, consisting of gravels and sands of the continental environment identified from the cuttings of the first tens to hundreds of metres for each well – this unit has been widely analysed in the literature (Regione Emilia Romagna and Eni-AGIP, 1998; Regione Lombardia and Eni-AGIP, 2002; Scardia et al., 2012; Ghielmi et al., 2013);

• late Pliocene–Pleistocene sand-rich sequences consisting of sands and clayey sands with clayey interlayers of deep-marine to continental environments (Sabbie di Asti Group), connected to the latest filling of the Po Plain foredeep;

• late Miocene–late Pliocene clastic deposits connected to the evolution of the northern Apennine foredeep and top-thrust basins, including

  • i.

    clay-rich units made up of variably silty clays with minor sand (Argille del Santerno Fm.)

  • ii.

    sand-rich (mostly foredeep turbiditic) units made up of thick sand beds with minor clay and thin-bedded sand–clay repetitions (e.g. Marnoso-Arenacea Fm., Bagnolo Fm., Fusignano Fm., Caviaga Fm., Porto Garibaldi Fm., Porto Corsini Fm.)

  • iii.

    conglomeratic units made up of shallow-marine and fluvio-deltaic conglomerates and sands with minor clay (e.g. Sergnano Fm., Cortemaggiore Fm., Boreca Conglomerate);

• early–late Miocene marly sequences, consisting of marl, clayey marl, clay, and sandy marl with sand intercalations recording deposition on the foreland; ramp (e.g. Marne di Gallare Group)

• Triassic to Eocene undifferentiated carbonate units consisting of prevailing limestone (mainly mudstone and packstone–grainstone) and dolomite with subordinate marl;

• Variscan crystalline basement.

We collected additional data related to the lithological characteristics of the subsurface sedimentary units from cutting descriptions. The cutting descriptions contain the information collected during mud logging, where rock fragments from the borehole reach the surface due to the circulation of the drilling fluids. Those data were combined with the SP and GR logs and with lithological data from core sample analysis reported in the well profile to characterize the lithology of the entire well. In total, we identified nine macro-lithologies, which are listed below together with the descriptions commonly found in the well profiles:

• gravel (e.g. polygenic gravel, gravel prevalent, polygenic gravel and sand with shale interbeds, and gravel and pebbles with sand interbeds)

• sand (e.g. sand, sand prevalent, shaly sand, fine sand, sand with shale interbeds, sand and shaly sands, and sand banks)

• cemented sand (e.g. cemented sand, fine-grained cemented sand, cemented sand and pebbles, calcareous cemented sand, and sand of variable cementation)

• shale (e.g. shale, silty shale, grey shale, marl shale, and prevalent shale)

• sand and shale alternances (e.g. grey shales and sand and shales with sand interbeds)

• conglomerates (e.g. polygenic conglomerate, polygenic conglomerate with shale, and polygenic conglomerate with sand interbeds)

• marl (e.g. marls, silty marls, grey silty marls, marls and sandy marls, and grey marls with cemented sand)

• dolomite (e.g. dolomite, calcareous dolomite, crystalline dolomite, shaly dolomite, and grey dolomite)

• limestone (mudstone, wackestone, packstone, and grainstone).

Well data represent the main constraint for the subsequent 3D geological modelling phase. As shown in Fig. 4, the distribution of the wells in the area is not homogeneous, and some regions are characterized by low well density. For this reason, the collection of the other primitive data (i.e. geological cross-sections and maps) was focused in the area indicated by the white polygon in Fig. 4, excluding the most isolated wells.

We performed a preliminary analysis of the sonic velocity variations with depth using the collected data and the newly created lithological and mineralization logs (Fig. 8). The sonic logs display transit time values in the range of 40–140 µs ft⁻¹ (131–459 µs m⁻¹). However, we observed higher values (approx. 150–200 µs ft⁻¹ (492–656 µs m⁻¹)) concentrated to the first tens to a few hundred metres from the surface which are mostly connected to gravel lithologies characterized by poor consolidation (Fig. 8). In some cases, unusually high or low sonic log values with respect to the general data trend for a certain lithology can be ascribed to the presence of thin layers with different lithological characteristics, to possible borehole damage, or to the presence of fractures. The currently available information does not allow a clear interpretation of their causes and even their removal. Most of the lithologies show a gradual decrease in transit time with increasing depth (Fig. 8b–d), and at about 4 km depth, the transit time flattens out, showing rather constant values for higher depths, independently from the lithology. The continuous decrease of the P wave transit time with depth reflects the increasing compaction of the sediments due to the overburden weight. The limestone is an exception, showing a relatively constant transit time independently of the depth. The presence of gas significantly increases the transit time with respect to water-saturated rocks: the so-called gas effect is outlined as a rock density reduction (Fig. 8a).

The abundance of each lithology within the main units defined in this study significantly affects the average mechanical rock properties of the entire unit (Fig. S2 in the Supplement). One of the advantages of the analysis performed is the possibility of assigning specific mechanical properties for each lithology and for different depths. Knowing the amount of each lithology present in the geological units, we can achieve a more precise characterization of the subsurface geological layers, which is fundamental for improving the quality of the geomechanical simulations. Applications of the sonic log analysis in defining mechanical rock properties of the subsurface of the Po Plain area can be found in Benetatos et al. (2023a, b). In the first work, the authors propose a workflow for geomechanical simulations through seabed monitoring. A significant step in the workflow involves the rock mechanical characterization that is performed through the analysis of the sonic log data that are converted to dynamic Young's modulus values that are specific to the different lithologies and to various depths. The same approach is also followed by Benetatos et al. (2023b) to calculate the mechanical properties of the Argille del Santerno Fm. They used well log data and compared them to those derived from laboratory analysis of core samples. This comparative analysis revealed relations between well log and laboratory-derived properties and contributed to the understanding of the deformation behaviour of this important geological formation that extends across northern Italy.

4.2 Geological cross-sections

We collected 6341 km of published geological cross-sections (Cassano et al., 1986; Lindquist, 1999; Casero, 2004; Picotti et al., 2007; Fantoni and Franciosi, 2009; Toscani et al., 2009; Wilson et al., 2009; Boccaletti et al., 2011; Pola et al., 2014; ISPRA, 2015; Maesano et al., 2015; Turrini et al., 2015; Livani et al., 2018; Table 1). Several procedures were implemented for uploading, calibrating, and revising the geological cross-sections in a 3D environment and, finally, for digitizing the selected horizons. The location maps (i.e. maps displaying the section traces) and the cross-sections were graphically re-arranged and improved to reduce imperfections and errors due to the low quality and/or image distortion. We then georeferenced the location maps in the QGis environment, digitizing the relative section traces (Fig. 4). Based on the geographically oriented section traces, the cross-section images (raster format) were properly uploaded in a 3D environment. The geological cross-sections composed of segments with different orientations were cut into several parts and separately imported. Where necessary (i.e. location inconsistency), the geological cross-sections were repositioned using the intersections with other data such as wells, surface geology (e.g. geological boundaries, faults), other geological cross-sections, and orographic and hydrographic features (e.g. rivers).

The revised and georeferenced cross-sections were finally uploaded into a 3D project in the Petrel^® software. We digitized four geological horizons that roughly correspond to the boundaries of the main lithological discontinuities recognized through the well data analysis. The horizons are, from the oldest to the youngest, as follows: the top of the carbonate succession (more marly in its upper portion), which divides a mainly carbonate succession (below) by a mainly siliciclastic succession (above), and it does not represent a chronostratigraphic boundary; the Pliocene base; the Calabrian base; and the base of recent continental deposits. In the regions characterized by a vertical duplication of the same unit due to the effect of thrusting, we digitized the hanging wall of the units until the hanging wall cut-off and then passed on directly to the footwall of the same unit. In such a way, a marked artificial step was generated; however, we deem this approximation to be necessary for the successive 3D modelling dataset creation, which does not integrate any fault element. In the data collection process, some public geological cross-sections were excluded from the database and the digitization process. For instance, the AGIP geological cross-sections reported in the Italian geological maps at 1:100 000 scale (sheets 75, 77, and 88; https://www.isprambiente.gov.it/it/attivita/suolo-e-territorio/cartografia/carte-geologiche-e-geotematiche/carta-geologica-alla-scala-1-a-100000, last access: 11 May 2021) were not used due to the geological interpretation based on scarce deep-data information (many deep wells and seismic reflection profiles were unavailable at that time), strongly differing from interpretations recently proposed on the basis of a large amount of subsoil data. Furthermore, some vintage sections (i.e. Bally et al., 1986) were excluded since some more recent geological cross-sections passing close to their traces show a more precise localization and better graphical properties, allowing for a more accurate digitization process.

For each digitized horizon, a delimited text file in ASCII format reporting the xyz coordinates is generated. Since the digitization process was performed manually, the data are provided with an irregular sampling step.

4.3 Geological maps

The 3D database also includes data derived from 10 published subsurface geological maps. We digitized the following maps:

• two isobath maps of the base of the recent continental deposits (QC1 horizon map by ISPRA, 2015; Scardia et al., 2012)

• one isobath map of the Calabrian base (QM1 horizon map by ISPRA, 2015)

• two isobath maps of the Pliocene base (Bigi et al., 1990a, b; ISPRA, 2015)

• four isobath maps of the carbonate succession top (Casero et al., 1990; Nicolich et al., 2004; ISPRA, 2015; D'Ambrogi et al., 2023)

• one isobath map of the magnetic basement top (Cassano et al., 1986).

The maps, available in hard copy and/or in raster format, were graphically revised and re-arranged in order to reduce possible scanning defects and distortions. The maps were then georeferenced, and the contour lines were digitized using QGis software. The digitized contour lines were exported and uploaded into the 3D Petrel^® database. We then verified in a 3D environment their consistency with other subsurface geological information, especially well data.

For each map, a delimited text file in ASCII format reporting the xyz coordinates of the contour lines was generated. Since the digitization process was performed manually, the data are provided with an irregular sampling step.

5 Data accuracy

A widespread accepted standard approach to address map accuracy is still not available. Recent studies suggest methods and standards for error analysis of geological or subsurface maps (e.g. Kint et al., 2020, and references therein). However, these methods are still affected by significant bias as they depend on the data availability and the criteria of data selection. Quality flagging is the basic approach used to quantify uncertainty within a spatial dataset and is done by assessing metadata fields. This method can be limited to indicate the presence or absence of data, or it can be very complex, producing a full range of quantitative error ranges (e.g. Bardossy and Fodor, 2001). Kint et al. (2020) presented an approach to assess data uncertainty for a well dataset in the Quaternary succession of the Belgian Continental Shelf. They produced confidence maps based on datasets from different origins and time periods. Their method consists of the following: (i) determination of the data density (how much data contribute to each grid cell to provide information on lateral and depth variability), (ii) direct mapping of measured errors and accuracies, (iii) transformation of the measured values or categorical quality flags into uncertainty percentages, and (iv) selection of data subsets based on the uncertainty maps. Not all these points are always feasible or necessary. This renders the method non-general. For example, in the case of few data or datasets without associated uncertainty, steps (ii) and (iv) are not recommended or feasible. Furthermore, the uncertainty drags errors due to the fact that instrumental (absolute accuracy, positioning accuracy) and human (expert judgement, data selection, data origin and representation) accuracy are often not quantifiable.

The only universal and dataset-unrelated rule when considering the geographic space comes from the first law of geography (Tobler, 1970) that states “everything is related to everything else, but near things are more related than distant things”. This law is the foundation of the concepts of spatial dependence and spatial auto-correlation and is used specifically for the inverse distance weighting method for spatial interpolation (Shepard, 1968).

Generalizing the approach proposed by Kint et al. (2020) with an application to arbitrary spatial data and using the inverse distance weight (IDW) for our analysis, we implemented, in a preliminary way, a method to weight the accuracy and/or confidence of geological surfaces. For each horizon, we quantified the following: (i) the data density contributing to the assessment of the lateral accuracy and the depth variability; (ii) the accuracy based on the data density and the spatial autocorrelation converted into a probability (inverse distance weight – IDW), describing the confidence of the data at each point of the study area; and (iii) the error associated with the depth of the geological surfaces due to discrepancies between the data of different origins where different guesses exist. Hence, for each point of the study area, we provided a value indicating the accuracy and a value indicating the estimated error of the depth of the geological surface.

Figure 9 (a) Grid nodes used to evaluate the data accuracy and uncertainties for geological surfaces presented for the study area. At each node, based on spatial autocorrelation, we calculated the inverse distance weight (IDW) with respect to point observations (checkpoints). (b) The model adopted for the IDW describes how uncertainties increase, moving away from the best guess at checkpoints (i.e. at wells, along geologic cross-sections, or on depth maps); at the same distance from a checkpoint, well data give higher IDWs (red points) in comparison to data derived from geological cross-sections or interpolated maps (green and black points).

Table 5 Parameters of each unit. Checkpoint refers to the number of locations where the data are observed. P is the checkpoint order, where p=1 for well data, p=2 for sections, and p=3 for maps.

Figure 10 IDW values at grid nodes calculated for each surface considering distances from well data only (dashed lines) or integrating distances from well data plus geological cross-section plus geological maps (solid lines). The area between the two lines for each proposed surface is larger for greater uncertainties. Note that the uncertainties increase with the depth of the surfaces and that the total number of checkpoints does not always increase the confidence level of the data (as, for example, in the recent continental deposits).

Figure 11 The spatial distribution of the IDW for the four analysed surfaces. The IDW values represent the confidence of the variable of interest (level depth) at each node of the grid (1 indicates max confidence).

To apply our model, we started by edging the study area including the observation region and covering it through an evaluation grid (Fig. 9a). In surface and subsurface mapping, the observations come from properly identified sites. The locations of these sites (hereafter checkpoints) are known and tag the observations. The variable of interest (i.e. horizon depth) exists in every point of the region (i.e. grid nodes) but is observed only in a finite set of locations (checkpoints). The variability model describes how uncertainties increase moving away from the checkpoint with respect to the best guess, built according to the observations and the geological constraints. In our case study, the properly identified sites are points (wells), lines (the traces of geological cross-sections digitized at discrete points), and polygons (here intended as raster fields digitized at discrete points, i.e. isobaths of geological maps). Note that not all kinds of checkpoints can be considered at the same level of confidence since well data give better depth constraints to those derived from geological cross-sections or interpolated maps. For this reason, it is convenient to assign higher specific weights to first-order checkpoints (well data in our case) with respect to higher-order checkpoints (sections or maps). Once the hierarchical subdivision has been made, the inverse distance (ID) principle is used to model the uncertainties at each point of the study area (i.e. at grid nodes) based on spatial autocorrelation with respect to checkpoints:

$$\begin{array}{}\text{(1)}& \text{ID}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+r^p}},\end{array}$$

where ID is the inverse distance, r is the distance between the point of observation and the checkpoint, and p is the checkpoint order; in this study, p=1 for well data, p=2 for sections, and p=3 for maps. The inverse distance (ID) is then normalized to obtain the inverse distance weight IDW as follows:

$$\begin{array}{}\text{(2)}& \text{IDW}={\displaystyle \frac{\text{ID}-{\text{ID}}_{\min }}{{\text{ID}}_{\max }-{\text{ID}}_{\min }}}.\end{array}$$

The IDW assigns a confidence based on distance autocorrelation that decreases more gently when considering first-order checkpoints (red curve goes under 0.2 at a distance of 7 km; Fig. 9b) with respect to higher-order checkpoints (green and black curves go under 0.2 at a distance of 2 km; Fig. 9b). The total weight calculated at each point of the study area is the mean value of the weights calculated with respect to different-order checkpoints. This means that IDW = 1 is the best guess assigned to grid nodes that, in the ideal case, lie – at the same time – over checkpoints of order 1 to n (i.e. in our case, over a well crossed by a section and coinciding with an interpolation point of the map). Statistics of the dataset are reported in Table 5.

5.1 Data density

Four out of the five processed geological surfaces (i.e. the base of recent continental deposits, the Calabrian base, the Pliocene base, and the carbonate succession top) were reconstructed using checkpoints of order 1 (i.e. well data), 2 (i.e. geological cross-sections), and 3 (i.e. subsurface geological maps). The reconstruction of the magnetic basement top derives from a unique source (Cassano et al., 1986), and, for this reason, we avoided including this surface in the data accuracy analysis.

The analysis of our dataset (Table 5) shows that the three shallowest surfaces (i.e. recent continental deposits, Pleistocene base, and Pliocene base) are the most constrained with up to 139 checkpoints of the first order, whereas the carbonate succession top is constrained by a few checkpoints of order 1. The carbonate succession top and the Pliocene base are described by the highest amount of data, with a density larger than 74 total checkpoints to each grid cell.

5.2 Data analysis

The total number of checkpoints does not always increase the confidence level of the data. In fact, for the Calabrian base, the Pliocene base, and the carbonate succession top, the IDW values also increase by up to double the amount when calculated considering checkpoints of all orders, whereas they decrease for recent continental deposits (dashed lines in Fig. 10). In the latter case, higher confidence is obtained when considering only the order-1 checkpoints (solid lines in Fig. 10). The best-constrained surface is confirmed to be the Pliocene base, with a high number of both well data (i.e. order-1 checkpoints) and total data. The spatial distribution of the IDW for the proposed surfaces is reported in Fig. 11 and represents the confidence in the variable of interest (level depth) at each node of the grid (1 indicates maximum confidence). All the maps show maximum IDW ≈ 0.7 and a mean value of 0.06 (recent continental deposits – Fig. 11a), 0.09 (Calabrian base – Fig. 11b), 0.14 (Pliocene base – Fig. 11c), and 0.08 (carbonate top – Fig. 11d). The IDW values indicate that the Pliocene surface has a confidence level almost double than that obtained for the other surfaces on average. Interactive figures (HTML format to be opened in a web browser) of each map reported in Fig. 11 are available in the Supplement (S1.zip) and provide detailed statistical information at each node of the interpolation grid.

5.3 Uncertainty associated with depth

The four surfaces analysed in this work derive from different origins and types of data (see Sect. 3; only the magnetic basement top derives from a unique source). The variability in quality, periods of creation, owners, and compiler sensibility (human error) of the datasets produces large – not quantifiable – uncertainty affecting the input data. Furthermore, since most of the collected public data derive from interpretations of seismic profiles, the data in depth are strongly influenced by the velocity model used for the depth conversion. Unfortunately, most of the primary data are not provided along with the used velocity model, and for this reason, it was not possible to take into account this variable. Further, the study area is non-uniformly described by the dataset.

Figure 12 Range of depth variations interpolated at each point of the evaluation grid based on the different datasets available for the four surfaces.

Figure 13 Geographic distribution of the range of depth variations for the four analysed surfaces. It is useful to compare these data with the IDW values in Fig. 11 in order to obtain more complete and comprehensive information.

Each grid node is associated with one to four values of depth depending on the number of available identified sites (i.e. checkpoints from different sources) coinciding with that node. Hence, it turns out that it is possible to calculate the maximum and minimum depth of each surface at each node if at least two depth values are available in that specific node. The range of depths at the nodes represents the uncertainty of the level description. Figure 12 shows the depth variation for the four surfaces calculated considering the maximum and minimum values for each node with at least two values. The bottoms of recent continental deposits (light-blue band) and of the Pleistocene unit (yellow band) show a small variation in depth (<3 km) in comparison to the Pliocene (orange band – up to 6 km) and carbonate (green band – up to 8 km) bases. The maps in Fig. 13 illustrate the geographic distribution of depth evaluation uncertainty. Interactive figures of each map reported in Fig. 13 are available in the Supplement (S1.zip).

6 Derived data: methods and results

The construction of an accurate 3D geological model requires good coverage and consistency of the source data. Improving the multi-source data consistency helps to avoid errors and distortions during the 3D model generation. To ensure the most significant data coverage and consistency in the area, we revised, compared, and finally integrated the collected primitive data (i.e. wells, geological cross-sections, and maps).

Based on the data accuracy analysis and the uncertainty associated with level depth (see Sect. 5), we totally or partially removed the digitized primitive data showing the highest discrepancy with other data. For this purpose, as mentioned in the previous section, we assigned a priority order to the data according to their type. The highest priority is attributed to well data as they give the best depth constraint of the units (order-1 checkpoints). Hence, the well data represent the main control points of our model, and none of them were excluded during the integration process. Regarding geological cross-sections and maps, if a difference in the depth of the surface is observed (Fig. 13), we kept those data showing the best fit with the well data. In other cases, we kept the most recent data, except when they were very discordant from the average data (for example, geological sections with a horizon depth very different from that indicated by other more reliable data) or from the predominant interpretative schemes (for example, geological cross-sections with a tectonic style completely different from the preponderant interpretations). In particular, it should be pointed out that the base of recent continental deposits is defined from well data and the isobath map reported in Scardia et al. (2012), which integrates the aquifer data of the Lombardy (Regione Lombardia and Eni-AGIP, 2002) and Emilia Romagna regions (Regione Emilia Romagna and Eni-AGIP, 1998). The inconsistent data removal was performed manually. These revised data integrated with the top of the magnetic basement, the only unit deriving from a single source (i.e. Cassano et al., 1986) and hence without any integration, were collected in ASCII format, reporting the xyz coordinates of the digitized data.

Figure 14 (a–e) Geological maps representing the gridded surfaces generated with the revised primitive data: (a) base of the recent continental deposits, (b) Calabrian base, (c) Pliocene base, (d) carbonate succession top, (e) magnetic basement top. (f) Example of the data digitized and the geological model in 3D view cut along a SW–NE-oriented section.

We eventually constructed a 3D geological model of the rock volume interposed between the magnetic basement top and the topographic surface (middle Triassic–present) by means of the revised primitive data. To define the top of the modelled rock volume, since our study area is located both offshore and onshore (Fig. 2), we joined the land topography, deriving from a public digital elevation model of the whole Italian territory at an original resolution of 10 m cell size (Tarquini et al., 2007), with the bathymetry of the offshore area (https://www.gebco.net/, last access: 15 March 2018). The 3D geological model is made up of unfaulted surfaces obtained by gridding the revised and integrated multi-source primitive data. For the gridding process, we applied in Petrel^® software a convergent interpolation algorithm with a 250 m sampling step. After gridding, we quality-checked and, where necessary, manually improved the modelled surfaces. The 3D model does not integrate any fault element. Thus, the modelled surfaces do not show any dislocation but rather shows sudden height differences in correspondence with the major tectonic structures, mostly where thrust sheets produce the vertical superimposition of the same unit (Fig. 14). The accuracy of the derived geological surfaces mainly depends on the quality and the quantity of source data (i.e. primitive data).

The aforementioned gridded surfaces represent the main boundaries that define units with different mechanical behaviour and subdivide the model into the following five sub-volumes, from top to bottom: the continental portion of the Quaternary deposits, the Pleistocene Asti sands, the Pliocene sands (i.e. Santerno, Porto Corsini, and Porto Garibaldi sandstones), the upper Eocene to Messinian siliciclastic formations (i.e. Gonfolite-Gallare Marls, Marnoso Arenacea Fm., San Donà Marls, Glauconie di Cavanella Fm., Fusignano Fm., Sergnano gravels, etc.), and the carbonate-marly succession.

The base of the recent continental deposits (Fig. 14a) was reconstructed in a restricted area, where it was intercepted by wells and modelled in the literature (Scardia et al., 2012). Areas with recent prodelta, delta, and marine sediments have been excluded (Fig. 14a). The base of the recent continental deposits consists of a slightly articulated surface that generally deepens from the peripheral sectors of the model area to the median areas of the Po valley. It is characterized by a depocenter area located between the Emilia and Ferrara arcs, where it reaches its maximum depth in the Parma area. The accuracy of this surface is generally intermediate, with a strong constraint given by the isobath map digitized from Scardia et al. (2012) (Fig. 4).

The Calabrian base (Fig. 14b) generally delimits the lower contact between the Asti sands and the Gelasian and Pliocene sandstones, except where, due to the erosion affecting the anticline culminations (in correspondence with the Emilia and Ferrara arcs), the Pleistocene deposits are unconformably in contact with the Miocene ones (on the top of the Emilia Arc anticlines). This surface progressively deepens from the model borders to the inner part of the Po Plain, except in the Ferrara Arc area, where some culminations can be observed. It reaches its maximum depth in the Parma–Mantova area and in the coastal and northern Adriatic areas. The accuracy of this surface is variable but generally well constrained by most of the available geological cross-sections and numerous wells. Moreover, in the Brescia–Mantova area, it is constrained by an accurate depth contour map (ISPRA, 2015; Fig. 4c).

The Pliocene base covers almost the entire area except for the Emilia and Ferrara anticlines, where it was eroded due to the tectonic uplift (Fig. 14c). This lithological boundary shows a rather articulated morphology characterized by pronounced culmination (located on the main anticline axes of the Emilia and the Ferrara arcs) and depocenter areas (between the Emilia and Ferrara arcs, immediately south of the Ferrara Arc, and in the coastal area). The accuracy of this surface is variable, but it is generally constrained by a large number of wells, geological cross-sections, and two subsurface geological maps covering the entire model area (Bigi et al., 1990a, b) and the GEOMOL project study area (ISPRA, 2015; Fig. 4c).

The carbonate succession top separates the siliciclastic formations (upper Eocene–present) from the carbonate and marly ones (Triassic to middle–upper Eocene in age). This surface deepens from NE to SW (beneath the northern Apennine front), with a local and pronounced elevation at the Ferrara Arc, where carbonates were markedly uplifted by the Apennine orogenic process (Fig. 14d). The accuracy of this surface is variable. It is constrained by several geological cross-sections but, due to its considerable depth, only a few deep well data items, mostly located on the anticline culminations. Three geological maps located in the eastern portion of the model (Casero et al., 1990; Nicolich et al., 2004) and in the central portion of the model (D'Ambrogi et al., 2023) constrain the surface of the carbonate succession top (Fig. 4c).

As mentioned above, the gridded magnetic basement top (Fig. 14e) is based on the published depth contour map by Cassano et al. (1986); consequently, its accurateness depends only on the source data accuracy.

7 Data availability

Primitive and derived data are available for open access in ASCII format at the following link: https://doi.org/10.5281/zenodo.8126519 (Livani et al., 2023). The format and the organization of data are explained in the related description file.

The dataset can be imported into any software handling 3D geological data, well data information, and spatial vector data formats.

8 Conclusions

The database provides a collection of geological and geophysical data of the Po Plain subsurface and of the northern Adriatic Sea; these were collected and digitized from the literature and from open repositories. We digitized data from the composite logs of 160 wells drilled in the area. Borehole information (i.e. wellhead coordinates, rotary table elevation, measured depth, true depth, total depth, and deviation survey) and spontaneous potential, gamma ray, and sonic logs are provided in ASCII format. Five horizons were digitized from 61 geological cross-sections and from 10 geological maps, ordered as follows from the oldest to the youngest: the top of the magnetic basement, the top of the carbonate succession, the base of the Pliocene, the base of the early Pleistocene (i.e. near top Gelasian; see Sect. 3), and the base of recent continental deposits. The digitized data are provided in ASCII format, reporting the xyz coordinates of the digitized surfaces.

Through an accuracy analysis performed on the primitive data and their subsequent processing, a new set of data has been created (i.e. derived data). Since the primary data show a depth uncertainty, we accurately revised the primary data by integrating only the data showing the best fitting. From these data, we generated a simplified 3D geological model characterized by several gridded surfaces of the main geological units (Fig. 14f).

Through the elaboration of the digitized logs, it is possible to directly extract geophysical and mechanical properties of the rock volume interposed between the gridded surfaces (e.g. P wave velocity) and to obtain further derived properties. For example, sonic velocities can be converted to mechanical rock properties, such as Young's modulus or Poisson ratio, that find applicability in geomechanical simulations which are performed to evaluate the ground subsidence and/or uplift phenomena (Carminati and Di Donato, 1999; Benetatos et al., 2017; Antoncecchi et al., 2021) or the change of stress field in a specific area in response to natural processes or anthropic activities (e.g. Schutjens et al., 2010; Radwan and Sen, 2021; Hemami et al., 2021; Sangnimnuan et al., 2021).

The dataset developed in the present work supports applications in a wide range of research areas with benefits for scientists, practitioners, and decision makers. As an example, once populated with the values of seismic velocity, the 3D geological model can find several applications in seismological studies. It can be used to improve the procedure and to reduce the uncertainties during earthquake location, it can contribute to the calculation of more accurate focal mechanisms, and it can perform wave propagation and ground motion simulations (e.g. Magistrale et al., 1999; Süss et al., 2001; Molinari et al., 2015; Livani et al., 2022). The 3D model also represents a starting model in perturbation studies, such as linearized inversions of travel times for crustal velocities (e.g. Magistrale and Day, 1999) or for studies related to the seismic waveforms for crustal structure; moreover, it can be used to derive densities and to compare them to gravity observations (Roy and Clayton, 1999).

In conclusion, this database will be useful in better defining and mitigating the possible natural and anthropogenic risks for the preservation of the environment and the safeguarding of the social and economic interests of the territory, contributing to a better and more efficient management of subsoil resources.