3D Engineering Geological Modeling to Investigate a Liquefaction Site: An Example in Alluvial Holocene Sediments in the Po Plain, Italy

: Liquefaction-induced surface manifestations are the result of a complex geological–geotechnical phenomenon, driven by several controlling factors. We propose a multidisciplinary methodological approach, involving engineering geologists, geomorphologists, sedimentologists, and geotechnical engineers, to build a 3D engineering geological model for liquefaction assessment studies. The study area is Cavezzo (Po Plain, Italy), which is a municipality hit by superﬁcial liquefaction manifestations during the Emilia seismic crisis of May–June 2012. The site is characterized by a Holocene alluvial sequence of the ﬂoodplain, ﬂuvial channel, and crevasse splay deposits prone to liquefaction. The integration of different geotechnical investigations, such as boreholes, CPTm, CPTu, and laboratory tests, allowed us to recognize potentially liqueﬁable lithological units, crucial for hazard assessment studies. The resulting 3D engineering geological model reveals a strict correlation of co-seismic surface manifestations with buried silty sands and sandy silts within the shallow 10 m in ﬂuvial channel setting, which is capped and laterally conﬁned by clayey and silty deposits.


Introduction
Liquefaction has been one of the major causes of damage to structures and infrastructures in numerous early and recent earthquakes (e.g., Charleston, SC, USA, Mw 7.  The study area is located on the southern limb of the buried Mirandola antiform [20,21]. The bedrock consists of marls and sands of the Pliocene and lower Pleistocene "Argille Azzurre" formation and middle Pleistocene "Sabbie di Imola" formation [22]. The seismic hazard of the study area is associated with the tectonic activity of the buried "Ferrara folds". The recent tectonic activity of the buried faults [23] accounts for several drainage anomalies of the alluvial plain [24]. Cavezzo is located less than 10 km far from the epicentre of the 29 May 2012 earthquake (Mw 5.8) [25]. The shallow lithostratigraphic succession was deposited mostly by the Secchia River, whereas the deeper sands from the Po River [26].
During the 29 May 2012 earthquake [27,28], co-seismic effects, such as liquefaction, water level rise, and fractures, have been observed along the buried abandoned channel of the Secchia River (Figure 1b), are located along this riverbed, which was active during Roman and Medieval times, until the 13th century.

Hydrogeological Setting
The areas most vulnerable to earthquake-induced liquefaction phenomena, according to historical information, are fully saturated, loose sandy deposits [29,30]. From the hydrogeological point of view, the subsoil of Cavezzo is mainly constituted of aquitard or aquiclude silty-clayey sequences, with layers of sand [31]. The water level of the shallow phreatic aquifer in Cavezzo is continuously monitored by three phreatimeters (Figure 2; codes , located in the first 3 m from the ground level, acquiring measurements every eight days, and managed by ARPAE Emilia Romagna (http://cloud.consorziocer.it/FaldaNET/retefalda/index, accessed on 19 December 2021). The study area is located on the southern limb of the buried Mirandola antiform [20,21]. The bedrock consists of marls and sands of the Pliocene and lower Pleistocene "Argille Azzurre" formation and middle Pleistocene "Sabbie di Imola" formation [22]. The seismic hazard of the study area is associated with the tectonic activity of the buried "Ferrara folds". The recent tectonic activity of the buried faults [23] accounts for several drainage anomalies of the alluvial plain [24]. Cavezzo is located less than 10 km far from the epicentre of the 29 May 2012 earthquake (Mw 5.8) [25]. The shallow lithostratigraphic succession was deposited mostly by the Secchia River, whereas the deeper sands from the Po River [26].
During the 29 May 2012 earthquake [27,28], co-seismic effects, such as liquefaction, water level rise, and fractures, have been observed along the buried abandoned channel of the Secchia River (Figure 1b), are located along this riverbed, which was active during Roman and Medieval times, until the 13th century.

Hydrogeological Setting
The areas most vulnerable to earthquake-induced liquefaction phenomena, according to historical information, are fully saturated, loose sandy deposits [29,30]. From the hydrogeological point of view, the subsoil of Cavezzo is mainly constituted of aquitard or aquiclude silty-clayey sequences, with layers of sand [31]. The water level of the shallow phreatic aquifer in Cavezzo is continuously monitored by three phreatimeters (Figure 2; codes , located in the first 3 m from the ground level, acquiring measurements every eight days, and managed by ARPAE Emilia Romagna (http: //cloud.consorziocer.it/FaldaNET/retefalda/index, accessed on 19 December 2021).  The visual inspection of the phreatic level trend shows annual cyclic water leve tuations of about 2.5 m, with a minimum in September and October and maxim March and April (Figure 2a). These cyclic fluctuations are related to rainfall, but the groundwater pumping for agricultural purposes [32]. Therefore, two field cam were planned to measure the water level depth in private wells. Considering that t riod of minimum and maximum depth corresponds to higher and lower risk for liq tion susceptibility, these campaigns were conducted in March and September The visual inspection of the phreatic level trend shows annual cyclic water level fluctuations of about 2.5 m, with a minimum in September and October and maximum in March and April (Figure 2a). These cyclic fluctuations are related to rainfall, but also to the groundwater pumping for agricultural purposes [32]. Therefore, two field campaigns were planned to measure the water level depth in private wells. Considering that the period of minimum and maximum depth corresponds to higher and lower risk for liquefaction susceptibility, these campaigns were conducted in March and September 2018, respectively. Then, the water level depth measurements (36 and 23, respectively), have been interpolated using the Kriging approach of Geostatistics [33]. The maps of the water level depth show lower values in the proximity of the Secchia River (1 m, in March 2018) and higher values in the northern part of the study area (4-5 m deep; Figure 2). These data indicate a higher liquefaction risk in the southern part of the study area.

Geotechnical Investigations
For the geotechnical characterization of the site, data from previous field campaigns were collected, and new investigations were performed [34]. Overall, the data consist of borehole logs, mechanic, electric, seismic cone penetration tests (CPTm CPTu, and SCPTu), and laboratory tests, such as grain size distribution curves and Atterberg limits. In particular, the collected data belong to different databases that are listed as follows:

Borehole and Grain-Size Distribution Analysis
We analysed boreholes 9 to 40 m deep: eight drilled along the Secchia River paleochannel, three within the artificial levee of the active riverbed, and one within the alluvial plain ( Figure 1). Figure 3 shows the location of the sampled sand boils and boreholes, analysed by the engineering and geotechnical laboratories of the University of Pavia and Modena. Additional information is provided in Table S1, in the Supplementary Materials.
Grain-size analyses of sand boils, ejected near the borehole 16-138-S1 (Uccivello sand boil, see the location in Figure 3), were compared with the sand from the subsurface ( Figure 5).

Cone Penetration Tests
The mechanical and electrical cone penetration test (CPTm and CPTu) results were used to draw the stratigraphic logs to build the 3D engineering geological model. The CPTu provides pore water pressure (u) values, sampled every 2 cm, whereas the CPTm is sampled every 20 cm.
Thus, the mechanical and electrical cone penetration test interpretation is performed using the classifications of Robertson (1990) [35] and Schmertmann (1978) [36] for CPTu and CPTm, respectively.         The chart of Robertson (1990) [35] allows for distinguishing nine soil types (normalized soil behaviour type, SBTn classes). The radius of concentric circles of the chart can be used as soil behaviour type index, Ic. The Schmertmann (1978) [36] approach allows us to classify of different lithotypes. Geologismiky Geotechnical, CPeT-IT v 1.6, and GeoStru-Static Probing software were used to analyse the CPTu and CPTm measurements, respectively.

Cone Penetration Tests
The mechanical and electrical cone penetration test (CPTm and CPTu) results were used to draw the stratigraphic logs to build the 3D engineering geological model. The CPTu provides pore water pressure (u) values, sampled every 2 cm, whereas the CPTm is sampled every 20 cm.
The chart of Robertson (1990) [35] allows for distinguishing nine soil types (normalized soil behaviour type, SBTn classes). The radius of concentric circles of the chart can be used as soil behaviour type index, Ic. The Schmertmann (1978) [36] approach allows us to classify of different lithotypes. Geologismiky Geotechnical, CPeT-IT v 1.6, and GeoStru-Static Probing software were used to analyse the CPTu and CPTm measurements, respectively.

LC from Boreholes Soil Classification from Schmertmann (1978) [36] Ic Range SBTn from Robertson (1990) [37]
Clay with peat (At) Organic clay and mixed soils Non-liquefiable clay with peat soils Ic > 3. 5  Among the seven dominant lithologies, the clayey silt (La) was included with the lithological class of the clay (A), being not liquefiable soils, and the lithological class of the clayey sandy silt (Las) was included in the sandy silt (Ls).
The anthropogenic deposits were distinguished using the boreholes and cone penetration tests. The thickness of these deposits was estimated, considering the depth interval with qc and fs is higher than the other values measured along with the entire profile, representing anomalous measurements. An example is reported in Figure 7.
Among the seven dominant lithologies, the clayey silt (La) was included with the lithological class of the clay (A), being not liquefiable soils, and the lithological class of the clayey sandy silt (Las) was included in the sandy silt (Ls).
The anthropogenic deposits were distinguished using the boreholes and cone penetration tests. The thickness of these deposits was estimated, considering the depth interval with qc and fs is higher than the other values measured along with the entire profile, representing anomalous measurements. An example is reported in Figure 7.

Methodology
The developed methodological approach is aimed to reconstruct the 3D engineering geological model, in order to identify the sandy layers susceptible to liquefaction phenomena, using boreholes, cone penetration tests, and classification tests ( Figure 8). The outputs of the 3D engineering geological model are used as the starting point to analyse the subsurface architecture, in order to distinguish homogeneous areas for liquefaction hazard assessment. The solid modeling approach has been adopted to construct the 3D engineering geological modelling architecture of Cavezzo [38].

Methodology
The developed methodological approach is aimed to reconstruct the 3D engineering geological model, in order to identify the sandy layers susceptible to liquefaction phenomena, using boreholes, cone penetration tests, and classification tests ( Figure 8). The outputs of the 3D engineering geological model are used as the starting point to analyse the subsurface architecture, in order to distinguish homogeneous areas for liquefaction hazard assessment. The solid modeling approach has been adopted to construct the 3D engineering geological modelling architecture of Cavezzo [38]. The software used to develop the geological model is Groundwater Modelling System (GMS) Aquaveo. The methodology consists of three phases, described as follows: • Phase 1. Conceptual model The software used to develop the geological model is Groundwater Modelling System (GMS) Aquaveo. The methodology consists of three phases, described as follows: • Phase 1. Conceptual model The first phase is aimed to collect and harmonize the data. The stratigraphic profiles of the 30 m deep sequence, obtained using boreholes logs, mechanical, and electrical cone piezometric tests, are compiled using the detected lithological classes. The database includes the geographic coordinates of each stratigraphic profile, as well as the depth and thickness of each layer.

Conceptual Model
All available geotechnical investigations, including the boreholes and cone penetration tests, were exploited to determine the soil stratigraphy and identify the soil type. The stratigraphic profiles were used to build the architecture of the subsoil model of the top 30 m sequence. For the borehole log interpretation, layers thinner than 40 cm were included in the lithological classes, located at the top or bottom, with similar soil behaviour type values. Therefore, the resulting vertical resolution of the model is approximately 40 cm. The analysis of cone penetration test, performed in a 20 m zone from the boreholes, confirms that these deposits have soil behaviour type index values comparable with the sandy silt (Ls).
The comparisons between pairs of boreholes and cone penetration tests within a distance of 20 m give insight into the difference in the stratigraphic profiles obtained using boreholes and cone penetration tests. Indeed, as previously reported, the sandy mixtures are not well-identifiable [41,42]. In some cases, the sandy silt and silty sand layer, detected using boreholes, are classified as clayey soils (Figure 9).
For this reason, a correction of the Ic was applied to increase the reliability of the CPTu and improve the stratigraphic profiles obtained by these different geotechnical investigations. sandy silt (Ls).
The comparisons between pairs of boreholes and cone penetration tests within a distance of 20 m give insight into the difference in the stratigraphic profiles obtained using boreholes and cone penetration tests. Indeed, as previously reported, the sandy mixtures are not well-identifiable [41,42]. In some cases, the sandy silt and silty sand layer, detected using boreholes, are classified as clayey soils (Figure 9). For this reason, a correction of the Ic was applied to increase the reliability of the CPTu and improve the stratigraphic profiles obtained by these different geotechnical investigations.
Even though numerous subsurface investigations previously carried out in Cavezzo for urban planning purposes and post-earthquake reconstruction, no comprehensive studies to harmonize and interpret the results are available. In this work, the 3D engineering geological model was constructed by laboratory and in situ tests. Overall, the lithostratigraphic architecture of the top 30 m of the sequence shows that the study area is characterised by significant vertical and horizontal variations. The resulting simplified Even though numerous subsurface investigations previously carried out in Cavezzo for urban planning purposes and post-earthquake reconstruction, no comprehensive studies to harmonize and interpret the results are available. In this work, the 3D engineering geological model was constructed by laboratory and in situ tests. Overall, the lithostratigraphic architecture of the top 30 m of the sequence shows that the study area is characterised by significant vertical and horizontal variations. The resulting simplified conceptual model has been used to interpret the engineering properties of soils. In particular, five stratigraphic units were identified (Figure 10 In the area of the NE-SW fluvial ridge in the urban zone of Cavezzo (Figure 1), the analysis of the cone resistance (qc) and sleeve friction (fs) shows a high lateral and vertical lithological variability, due to the facies transition between channel and levee. The values of qc are not uniform for the shallow 5-10 m from the ground level, moving from the northern and southern portion of the fluvial ridge, made up of sandy layers, ranging from 1.5 to 5 MPa. From 15 to 20 m below the ground level, the qc values of the deeper sands are higher, ranging from 2 to 15 MPa.
-Unit A; heterogeneous deposits, lithological classes clayey silt, and clayey sandy silt (La, Las), with interbedded thin silty sand (Sl) layers, corresponding to the recent alluvial plain; -Unit B; lithological classes sand (S), silty sand (Sl), and sandy silt (Ls), corresponding to the fluvial channel, and these deposits are the source of most of the visible liquefaction effects; -Unit C; clay (A) and clay with peat (At), corresponding to the palustrine depositional environment; -Unit D; clay (A) of the ancient alluvial plain; -Unit E; dense sands of the ancient fluvial channel. In the abandoned riverbed at the plain level (Figure 1), shallow sandy layers show qc values similar to the fluvial ridge. The thickness of these sediments is lower than that of the sandy soils localized in the fluvial ridge. Few cone penetration tests are available to analyse the lithological variabilities of the crevasse splays. These geomorphological contexts display interbedded thin layers of sandy silt and silty sand, with qc values not very different, with respect to the values measured in the fluvial ridge and abandoned riverbed at the plain level. In the alluvial plain, sandy layers are present in the first two meters from the ground level, with values of qc ranging from 0.6 to 4 MPa.

Construction of a 3D Engineering Geological Model
The 3D engineering geological model of Cavezzo was constructed with a spatial resolution of 20 m, using 10 borehole logs and 164 stratigraphic profiles, obtained using 111 mechanical cone penetration, 32 electrical cone penetration, and 21 seismic cone penetration tests (Figure 1). The resulting model shows 24 horizons (Figure 11a). The first one is represented by the anthropogenic deposits, distributed over the entire area, reaching a maximum thickness of 2 m in the urban area. The 3D north-to-south and west-to-east cross-sections (Figure 11b) show a higher heterogeneity in the fluvial ridge, with lenses of sandy silt and silty sand, whereas continuous horizons are present in the floodplain. In the southern sector of the fluvial ridge, where liquefaction phenomena were observed (Figures 1 and 3), four main stratigraphic units are present (Figure 11c). Underneath the anthropogenic deposits, a superficial clayey silt and sandy silt layer, 2 to 4 m-thick, with a low percentage of sand (from 2 up to 13%), is present. The layer below is constituted by sand, sandy silt, and silty sand, with a thickness of 5-7 m. The sand content of these deposits ranges from 50 to 80%, as reported in Section 3.1. The third stratigraphic unit is characterised by clay with peat and clay. In the deepest unit, clayey deposits are present. Overall, the model displays a superficial coverage of sandy silt, with a higher content of fine with a thickness, ranging from 2 to 15 m. The higher thickness of these shallow deposits is observed in the floodplain. Below these deposits, sandy layers are present. These deposits are not laterally extended and mainly localized in the fluvial ridges and abandoned riverbeds at the plain level. The thickness of these soils ranges from about 4 to 6 m. The deeper deposits are constituted by clay and clay with peat, with the last one representing a marker layer, useful for the stratigraphic correlations.
Over 15-20 m from the ground level, fine sandy soils are present, corresponding to the ancient fluvial channel deposits.
Furthermore, the reliability of the 3D modelling has been assessed by performing a statistical analysis of the spatial distribution of the geotechnical investigations. The kernel density analysis has been implemented to estimate the density of the geotechnical investigations in a neighborhood around these points [43]. A total of 33% of the study area is characterized by values from moderate to high density, and it is evident that the density of the data used to create the 3D engineering geological model is higher where the liquefaction phenomena were observed, and the geological setting is more complex, due to the deposits of the buried paleochannel of the Secchia River ( Figure 12). Overall, the model displays a superficial coverage of sandy silt, with a higher content of fine with a thickness, ranging from 2 to 15 m. The higher thickness of these shallow deposits is observed in the floodplain. Below these deposits, sandy layers are present. These deposits are not laterally extended and mainly localized in the fluvial ridges and abandoned riverbeds at the plain level. The thickness of these soils ranges from about 4 to 6 m. The deeper deposits are constituted by clay and clay with peat, with the last one representing a marker layer, useful for the stratigraphic correlations.
Over 15-20 m from the ground level, fine sandy soils are present, corresponding to the ancient fluvial channel deposits.
Furthermore, the reliability of the 3D modelling has been assessed by performing a statistical analysis of the spatial distribution of the geotechnical investigations. The kernel density analysis has been implemented to estimate the density of the geotechnical investigations in a neighborhood around these points [43]. A total of 33% of the study area is characterized by values from moderate to high density, and it is evident that the density of the data used to create the 3D engineering geological model is higher where the liquefaction phenomena were observed, and the geological setting is more complex, due to the deposits of the buried paleochannel of the Secchia River ( Figure 12).

Engineering Geological Units for Liquefaction Hazard Assessment (MOPS)
The analysis shows that the study area is characterized by nine MOPS. The liquefiable layers were identified in all the MOPS, except for MOPS 9, corresponding to the alluvial plain sector (Table 2). These liquefiable soils include sandy silts or silty sands, as well as sands mainly localized in the shallow 2-15 m from the ground level. Although the liquefiable layers are evident in all the MOPS, field evidence and observed co-seismic effects, after the 2012 earthquake, indicate that only MOPS 5, 6, and 7 were affected by liquefaction ( Figure 13). In MOPS 5 and 7, silty sands and sandy silts are located between 2 and 12 m from ground level, in most cases capped by a clayey soil of about 1 m thick. Clayey deposits are not superimposed on the liquefiable deposits (such as silty sands and sandy

Engineering Geological Units for Liquefaction Hazard Assessment (MOPS)
The analysis shows that the study area is characterized by nine MOPS. The liquefiable layers were identified in all the MOPS, except for MOPS 9, corresponding to the alluvial plain sector (Table 2). These liquefiable soils include sandy silts or silty sands, as well as sands mainly localized in the shallow 2-15 m from the ground level. Although the liquefiable layers are evident in all the MOPS, field evidence and observed co-seismic effects, after the 2012 earthquake, indicate that only MOPS 5, 6, and 7 were affected by liquefaction ( Figure 13). In MOPS 5 and 7, silty sands and sandy silts are located between 2 and 12 m from ground level, in most cases capped by a clayey soil of about 1 m thick. Clayey deposits are not superimposed on the liquefiable deposits (such as silty sands and sandy silts), in correspondence with the MOPS 6. Even if surface manifestations of liquefaction can be evident in areas with no low permeability cap present, the 3D engineering geological model of Cavezzo displays that, in this case, the field evidence of liquefaction are observed in correspondence of silty sands and sandy silts layers in MOPS 5, 6, and 7, which are confined by clayey and silty deposits, both vertically and laterally, forming bodies that were extended not more than a few hundreds of meters in the lateral direction. These results agree with previous findings for the co-seismic liquefaction phenomenon observed in other sites in Emilia Romagna [27,[44][45][46][47]. The sandy silts and silty sands layers, identified in the other MOPS, where no liquefaction phenomena were detected, are continuous bodies not confined by clayey-silty deposits.

Depth of the Water Level (m) from Ground Level (March 2018)
1 Liquefiable sandy silt layers, between 2 and 9 m from ground level.

Discussion
In this section, the contribution of the 3D engineering geological model for liquefaction studies is discussed. In particular, the 3D engineering geological model, reconstructed for Cavezzo, was used to support the identification of the liquefaction source layer and mechanisms that may control the ejection to the surface.

Influence of the Sedimentary Architecture on the Surface Manifestations of Liquefaction
Previous authors observed that the surface manifestations of earthquake-induced liquefaction occur when a base layer of liquefiable soils is overlaid by a confining layer, known as the cap or crust [6,47,48]. In particular, Ishihara (1985) [6] introduced an empirical approach, based on the relationship between the thickness of the liquefiable layer and the thickness of the overlapping non-liquefiable layer, to assess the potential for ground manifestation of liquefaction. The author analysed the liquefaction observed after three large earthquakes occurred in China and Japan, with Mw > 7.5. In this Section, the 3D engineering geological model of Cavezzo was exploited to extract the thickness of the liquefiable and non-liquefiable layer, in order to understand the mechanisms that may have controlled the surface manifestations of liquefaction. The superficial sandy silt, with a very low content of sand, is considered as the potential non-liquefiable layer, and the underneath silty sand is identified as the liquefiable layer. A grid of 20 m was created, and a centroid was carried out for each cell of the study area. For each point, the values of the thickness of the non-liquefiable (H1) and liquefiable (H2) layers were assigned. Layer thickness data were plotted in the graph of Ishihara (1985), also reporting the points of different types of observed co-seismic effects ( Figure 14). Most of the observed liquefaction and water level rise, as well as the fractures and liquefaction phenomena, are correctly predicted by the bounds introduced by Ishihara (1985), except for some liquefaction sites. These misclassified points could be due to the presence of laterally drifted conduits of the sand boils in the area with higher thickness of the liquefiable soils (Youd and Garris, 1995). The results are consistent for the area of the Uccivello sand boils localized in the southern zone of the fluvial ridge, aligned along the NE-SW direction. Therefore, the ground surface manifestations observed in the studied area seem to be related to sandy layers of about 1-3 m thick in the upper 10 m from the ground level, laterally and vertically confined.
ical approach, based on the relationship between the thickness of the liquefiable layer and the thickness of the overlapping non-liquefiable layer, to assess the potential for ground manifestation of liquefaction. The author analysed the liquefaction observed after three large earthquakes occurred in China and Japan, with Mw > 7.5. In this Section, the 3D engineering geological model of Cavezzo was exploited to extract the thickness of the liquefiable and non-liquefiable layer, in order to understand the mechanisms that may have controlled the surface manifestations of liquefaction. The superficial sandy silt, with a very low content of sand, is considered as the potential non-liquefiable layer, and the underneath silty sand is identified as the liquefiable layer. A grid of 20 m was created, and a centroid was carried out for each cell of the study area. For each point, the values of the thickness of the non-liquefiable (H1) and liquefiable (H2) layers were assigned. Layer thickness data were plotted in the graph of Ishihara (1985), also reporting the points of different types of observed co-seismic effects ( Figure 14). Most of the observed liquefaction and water level rise, as well as the fractures and liquefaction phenomena, are correctly predicted by the bounds introduced by Ishihara (1985), except for some liquefaction sites. These misclassified points could be due to the presence of laterally drifted conduits of the sand boils in the area with higher thickness of the liquefiable soils (Youd and Garris, 1995). The results are consistent for the area of the Uccivello sand boils localized in the southern zone of the fluvial ridge, aligned along the NE-SW direction. Therefore, the ground surface manifestations observed in the studied area seem to be related to sandy layers of about 1-3 m thick in the upper 10 m from the ground level, laterally and vertically confined. Figure 14. The thickness of liquefiable (H2) versus the thickness of the non-liquefiable (H1) layers for Cavezzo distinguishing the non-occurrence points and ground surface manifestations of liquefaction. The bounds introduced by Ishihara (1985) [6] are also reported.

Composition of the Liquefied Layers
The examined sands from boreholes and sand boils have an overall lithoarenitic composition, distinctive of the Secchia River [49]. The lithic association includes sedimentary fine-grained siliciclastic grains (siltstones and shales) and carbonate lithics (largely micritic limestones and calcite spars). Shales are well-lithified and -rounded, with an evident iso-orientation of the clay minerals; for these characteristics, they appear to have a detrital origin, derived from the erosion of the pelitic successions cropping out in the Northern Apennines. The comparison of the composition of the sand boil sand with the buried sands indicates a close similarity with the shallow sands, down to the depth of 4.5 m. These shallow sands are discriminated from the deeper ones, which are richer in carbonate grain content.

Conclusions
The paper proposes an innovative multidisciplinary approach to understand the earthquake-induced liquefaction manifestations at the ground surface by constructing a 3D detailed engineering geological model. The proposed methodology allows us to facilitate the joint use of different direct and indirect surveys, such as boreholes, CPTm, and CPTu. The methodology allows to detect potentially liquefiable engineering geological units, essential for hazard assessment (MOPS), using the spatial variation of the behavioural index Ic, obtained from the interpretation of CPTm and CPTu test results. Indeed, the methodology allows us to identify the liquefiable layers and their geometry in a complex geomorphological and lithotechnical context, located in an alluvial plain, characterized by fluvial depositional bodies prone to liquefaction. Confined sandy silts and silty sands bodies, susceptible to liquefaction, were recognized, especially in areas where the depth of the water level is shallower. The 3D engineering geological model highlights the subsoil architecture of these bodies in the subsoil, which matches with the surface manifestations of liquefaction that were observed after the seismic sequence in 2012.
The results show that the depositional geometry should be considered to identify the prone areas to co-seismic surface manifestations. These results may be significant for the study of liquefaction in other similar contexts.