Liquefied sites of the 2012 Emilia earthquake: a comprehensive database of the geological and geotechnical features (Quaternary alluvial Po plain, Italy)

This paper presents a comprehensive geological and geotechnical study of the whole area affected by liquefaction following the 2012 Emilia earthquakes, including all the available information from the field reconnaissance surveys, in situ tests, and laboratory analyses. The compilation was performed at 120 liquefied sites to verify and validate the reliability of liquefaction charts in alluvial sediments, and to assess liquefaction induced by the 2012 seismic sequence in the Emilia plain. The results reveal a wide range of grain sizes (from clean sands to sandy silts) and compositional characteristics (quartz-rich to litharenitic) in the 2012 ejecta, and show a strong relationship between the liquefaction and stratigraphic architecture of the subsurface. The availability of in situ tests at the liquefied sites makes it possible to verify and validate the reliability of the liquefaction charts in alluvial sediments with respect to the real observations. For the analyzed Emilia case studies, the use of non-liquefiable crust provides better estimations of the liquefaction manifestations when coupled with the thickness of the liquefiable layer rather than with the liquefaction potential index. Altogether, this work makes available to the international scientific community a consistent liquefaction database for in-depth earthquake studies.


Introduction
During the last decades liquefaction case studies have been reported worldwide as a result of earthquakes with estimated magnitude M ≥ 5.4 (Maurer et al. 2015), providing detailed information on liquefaction features and geotechnical properties (e.g. Suzuki et al. 2003;Cao et al. 2011;Green et al. 2014;Facciorusso et al. 2015;Wood et al. 2017). Data from recent earthquakes flowed into an open-source global database, the Next Generation Liquefaction project (Zimmaro et al. 2019). This project aims to improve the accessibility of coseismic effect archives and support the development of updated deterministic and probabilistic liquefaction triggering curves using the "simplified procedure", as originally introduced by Seed and Idriss (1971). These "simplified methods" for liquefaction susceptibility assessment include procedures proposed by Seed et al. (1985), Robertson and Wride (1998), Andrus and Stokoe (2000), Juang et al. (2002Juang et al. ( , 2006, Cetin et al. (2004), Moss (2003), Moss et al. (2006), Idriss and Boulanger (2008), Kayen et al. (2013), Boulanger and Idriss (2014), and Stewart et al. (2016). In this context, the 2012 Emilia (Italy) seismic sequence is of great interest because of the widespread liquefaction phenomena in Holocene alluvial sediments that have been well documented through detailed geological and geotechnical surveys (Emergeo Working Group 2013, Papathanasiou 2012, Regione Emilia-Romagna 2012). Numerous studies have documented the liquefaction characteristics at various key sites (e.g. Amoroso et al. 2017Amoroso et al. , 2020Lai et al. 2020;Tonni et al. 2015;Fontana et al. 2019;Meisina et al. 2019;Facciorusso et al. 2015;Civico et al. 2015;Rollins et al. 2021), but a comprehensive study of liquefaction over the whole area affected by the 2012 Emilia earthquakes, that includes both geological and geotechnical aspects, is still lacking.
The aim of the present study is to provide an unprecedented compilation of the entire set of available information, including post-earthquake surveys, in situ tests, and laboratory analyses for 120 liquefied sites. The implications of this research are twofold: (1) at a general scale, the study makes it possible to verify and validate the reliability of liquefaction triggering charts for alluvial sediments, using a significant dataset that includes a wide grain size (from clean sands to sandy silts) and compositional (quartz-rich to litharenitic) spectrum of sediments; and (2) at a regional scale, the work represents a pivotal contribution to the liquefaction assessment in the Emilia plain.

Geological and seismological setting
In 2012, northern Italy was hit by two mainshocks and several large aftershocks ( Fig. 1; Pondrelli et al. 2012), generated by various thrust structures. The strongest shock occurred on May 20 th (moment magnitude, M w = 6.1), with the epicenter near the town of Finale Emilia in the Emilia-Romagna Region. Many aftershocks followed in the same area, up to the local magnitude M L = 5.1. The second largest earthquake (M w = 5.9) took place on May 29 th to the west of the first mainshock and close to the town of Mirandola. The 2012 seismic events caused many human casualties as well as widespread damage to buildings and infrastructure (Cultrera et al. 2014).
The epicentral area of the seismic sequence, located south of the Po River (Modena, Ferrara, Bologna and Mantova provinces), was affected by the liquefaction phenomena (Fig. 1). The Emilia alluvial plain is crossed by several rivers (Fig. 1), flowing from the Apennines into the Po River (Secchia, Panaro), or directly reaching the Adriatic Sea (Reno). The alluvial plain is the surface expression of the foredeep basin of the Apennine chain, where the combination of fast subsidence and strong sediment input generated very thick Plio-Pleistocene successions (Ghielmi et al. 2010). The study area corresponds to the buried frontal portion of the compressive ramp and flat structures of the Apennines (Toscani et al. 2009;Martelli et al. 2017). The active faults are associated with moderate seismic activity (Michetti et al. 2012), characterized by shallow epicenters, well documented through the last centuries (Locati et al. 2011;Guidoboni et al. 2018). However, taking into account the social and economic regional context, the related expected seismic risk is rather high. The most significant historic earthquake, that impacted the Ferrara area in November 1570 (VIII MCS, estimated moment magnitude M w = 5.5), induced sand liquefaction and ground fracturing, both within the town of Ferrara (Guidoboni et al. 2018) and in the Reno River sediments at San Carlo, that also exhibited liquefaction in 2012 (Caputo et al. 2016).
Both the 2012 mainshocks triggered widespread liquefaction of granular sediments, sand boils, ground fractures, and gravitational lateral spreading. The liquefaction phenomena in the eastern sector of San Felice sul Panaro were mainly associated with the mainshock of the 2012 May 20 th , whereas those in the western sector were coupled with the May 29 th earthquake. It is worthy of note that at some sites liquefaction phenomena (e.g. San Felice sul Panaro) have been observed following both the May 20 th and 29 th shocks (Pizzi and Scisciani 2012;Emergeo Working Group 2013). Moreover, the May 20 th aftershocks, occurred within a close time interval (less than 4 minutes with M L between 4.8 and  Paolucci et al. (2015) for the M W = 5.9 on May 29 th 5.0), may have had some effects on pore water pressure build-up within the saturated cohesionless layers, especially in the municipality of Terre del Reno (villages of Sant'Agostino, San Carlo, Mirabello) where the most and largest liquefaction effects were observed (Sinatra and Foti 2015;Facciorusso et al. 2016).
Liquefaction mainly took place within fluvial channel deposits. A strong relationship exists between the siting of the liquefaction and the stratigraphic architecture of the subsurface. The liquefaction sites are typically confined to elongated strips, corresponding to fluvial buried sandy bodies of Holocene age (Civico et al. 2015). The sites analyzed in the present research consider the geotechnical investigations available approximately in the upper 30 m of depth, where the stratigraphic interval consists of late Pleistocene and Holocene deposits. The interval is divided into two superposed subsynthems: the upper Pleistocene AES7 and the Holocene AES8 (Stefani et al. 2018). The ejected liquefied sands sampled at the surface derive entirely from the Holocene units, in particular from sediments younger than roughly 4000 years.
The analyzed deposits record two main sediment inputs, fed by the Po River to the north, and by several streams flowing from the Apennine chain to the south (Secchia, Panaro and Reno rivers; Fig. 1). The Po River largely differs from the Apennines counterparts in its fluvial dynamics, grain size distribution, and petrographic composition. The two fluvial systems therefore generated sharply different sediments, in the northern and southern portions of the research areas, respectively.
During synglacial times (AES7 unit), the Po deposited large coarse sand bodies into braided river environments, whereas during the Holocene (AES8), it sedimented finer grained sands into meandering channels. The Holocene meander channels often reworked the upper Pleistocene deposits, generating thick continuous bodies of sands. The depositional geometry and geomorphic expression of some of the younger meander bodies are visible at the surface (Figs. 2b and c). During the last millennia, the Po started to generate hanging channels, less curved in plan, interspaced with interfluvial depressions, where finer grained sediments accumulated.
The Apennine rivers also provided large volumes of sediments, but finer grained compared to the Po deposits. During the last glaciation, great volumes of silt and argillaceous silt accumulated, framing thin bodies of fluvial sands (AES7). During the Holocene, the rivers flowing from the Apennine valleys produced large volumes of clayey or silty sediments, often rich in peat and formed flood plains and freshwater marshes, framing isolated fluvial channel bodies, mainly consisting of silty sands. Throughout the last millennia, the boundary between the deposits of the Apennine rivers and of the Po River progressively migrated northward. In medieval times, the Apennine rivers were unable to flow into the Po and formed large inland delta systems (Fig. 2a). The inland distributary channels of the Secchia, Panaro and Reno rivers prograded into large freshwater marshes and shallow lakes. A significant portion of the liquefied sands analyzed in the southern part of the study region accumulated in these inland delta distributary channels, which often maintain a strong topographic expression, as elongated ridges (Fig. 2a).

Field survey and sampling
The detailed 2012 field reconnaissance survey identified different coseismic geological features associated with liquefaction or fracture/liquefaction. These features include single sand volcanoes, scattered vents, coalescent flat cones, sand infilled water wells, fountains and manholes, elongated and aligned multiple sand volcanoes and sand flows from open fractures. Examples of these liquefaction features are reported in the photos in Fig. 3.
A total of 120 sites were analyzed and sampled in the alluvial plain area located between the Po, Reno and Secchia rivers in the Ferrara, Modena, Bologna and Mantova provinces   are available at https:// istit uto. ingv. it/ images/ colla ne-edito riali/ misce llanea/ misce llanea-2012/ misce llane a16. pdf ( Fig. 1). For the liquefaction sites we report, when available, structural and morphological data including: (1) morphology, diameter and thickness of the sand boils; (2) spacing, including length and strike of the sand boil alignments and associated fractures. The thickness of the extruded sand was up to about 40 cm, the maximum observed diameter of individual sand volcanoes was 10 m, and the coalescent sand volcanoes along fractures extended for a maximum length of 50 m (Emergeo Working Group 2013). Samples of the liquefied extruded sand were collected for sedimentological analysis. The samples were preserved in plastic bags and their weight was commonly between 200 and 500 grams. Moreover, qualitative surveys of liquefaction-induced land damage related to 2012 Emilia earthquakes were performed by means of a reappraisal of field and aerial identification and characterization of liquefaction phenomena. The land damage estimation was carried out following the qualitative scale proposed by van Ballegooy et al. (2014): • Category 1 -No observed ground cracking or ejected liquefied material; • Category 2 -Minor ground cracking but not observed ejected liquefied material; • Category 3 -No lateral spreading but minor to moderate quantities of ejected material; • Category 4 -No lateral spreading but large quantities of ejected material; • Category 5 -Moderate to major lateral spreading; ejected material often observed; • Category 6 -Severe lateral spreading; ejected material often observed.
This estimation is clearly relative to the 2012 damage and Categories 1 and 2 were not considered since the 120 analyzed sites are all liquefied sites with ejected material. For each site a detailed report sheet has been compiled (see electronic supplement of this paper).

Grain size and compositional analyses
A total of 120 sand boil samples were analyzed using standard techniques. To reduce sampling bias due to possible segregation of particle sizes during the sand boil formation, for each site we sampled the complete vertical sequence next to the ejection point; the deposit normally did not exceed 10 cm in thickness. For comparison, two additional samples (San Martino 2b and San Martino 6bis; see the electronic supplement of this paper) were also collected from the fine-grained distal portion of the ejected deposits. These samples are not plotted in the graph of the Figs. 4 and 5, presented in the next sections of the work, as they are not representative of the complete ejected deposit.
Mechanical sieving was performed for the sandy fraction along with hydrometer analysis or laser diffractometry analysis for the fine-grained sediments. Sand samples consisting of a few hundred grams were washed with dilute H 2 O 2 to remove organic matter and were air dried and mechanically sieved for grain size and compositional analyses.
Compositional analyses were carried out on sand samples by point counting (300 grains for each thin section) under transmitted-light microscopy. Analyses were performed on the 0.125-0.250 mm fraction, according to the Gazzi-Dickinson method designed to minimize the dependence of the analysis on the grain size (Zuffa 1985) and to support a comparison with previous compositional studies performed in the area (Lugli et al. 2004(Lugli et al. , 2007. The aim of the compositional study was to attribute the ejected sand to the specific buried fluvial channels of the alluvial system (Po, Secchia, Panaro and Reno rivers). Sorting was calculated according to Folk and Ward (1957) using the GRADISTAT software.

Geotechnical investigation
The presented work benefitted greatly from the large dataset of previous geological-geotechnical subsurface investigations, developed by the Emilia-Romagna Regional Administration (https:// ambie nte. regio ne. emilia-romag na. it/ it/ geolo gia/ carto grafia/ webgis-banch edati/ banca-dati-geogn ostica). The available dataset, compiled since the early 1990s, provides more than 65,000 logs, throughout the alluvial plain, including more than 5,000 piezocone penetration tests (about 1,000 are located within the epicentral region of the 2012 earthquake), deriving from microzonation studies and reconstruction projects developed after the seismic crisis (Regione Emilia-Romagna 2013; Martelli et al. 2013). The dataset of the 120 sites includes the available piezocone tests (CPTU), electrical cone penetration tests (CPTE), dilatometer tests (DMT), seismic dilatometer tests (SDMT) and boreholes with standard penetration tests (SPT) performed at the liquefied sites or in their proximity.

Data analysis
Relevant data for all 120 liquefaction sites has been assembled and analyzed. For each site, a comprehensive dataset has been developed including geographical and morphological data, liquefaction features, liquefaction-induced land damage and related seismic parameters. In addition, grain size and compositional analyses, paleoriver assignment, and geotechnical tests were assembled. Finally, liquefaction assessments have been performed using in situ tests. These data are organized as a single report sheet presented in the electronic supplement of this paper and are summarized in Table 1.

Grain size analysis
Sand blows show variable grain size distributions as illustrated in the cumulative curves of Fig. 4 and in Table 1. The index properties for each sand blow were estimated using the grain size distribution, and include (Table 1): • the fines content (FC) which represents the percentage of particles finer than 0.075 mm; • the coefficient of uniformity (U), given by the equation: where D x is the diameter of the generic xth percentile of the grain size curve; • the coefficient of curvature (C), given by the equation: • the sorting (σ), given by the equation: (2) C = (D 30 ) 2 D 10 ⋅ D 60        where x = Log 2 D x and D x is the diameter of the generic xth percentile of the grain size curve in millimeters.
The grain size curves and the calculated index properties were then used to determine the Unified Soil Classification System symbols (USCS classification, ASTM D2487-11 2011) reported in Table 1. Sand boils were most frequently classified as silty sand (SM) and sometimes as poorly graded sand with silt (SP-SM) or sandy silt (ML), considering that FC varies on average between 10 and 40%, U is typically below 30 and C varies between 1 and 4.
The coarser samples (from medium sand to fine sand) are from sand blows in the northern area, between Bondeno and Quistello (see Fig. 1 for the their location), that is dominated by the Po River paleochannels. The sand blows from the eastern sector, dominated by the Reno paleochannel between Sant'Agostino and Mirabello (see Fig. 1 for the their location), are generally finer grained compared to those of the northern sector, and show a larger spectrum from fine-grained sands to silty sands. The samples from the western area (San Felice sul Panaro, Concordia sulla Secchia and Cavezzo) are characterized by the largest percentage of fines (up to 50%) and consist of fine-grained silty sands and sandy silts. These sands were ejected in the area dominated by the Secchia and Panaro rivers.
The grain size ranges for the four main fluvial systems (Po, Reno, Secchia and Panaro) are outlined in Fig. 4, showing that the Secchia and Panaro river sands are generally finer compared to the Po sediments, while the Reno deposits cover a wider grain-size range. Figure 4 indicates also that the liquefied sediments are within the grain size range of deposits typically susceptible to liquefaction according to observations of Tsuchida and Hayashi (1971). As shown in Fig. 5a, sorting of ejected sand is moderate to poor for all of the source rivers and the lowest degree of sorting generally characterizes the samples with higher fines content (Fig. 5b).

Petrographic composition analyses
Sand blows and dikes have a variable composition ranging from (1) lithoarenite to (2) quartz-rich lithoarenite and (3) quartzarenite (Fig. 6). The lithoarenitic sands contain an abundant lithic association including sedimentary fine-grained siliciclastic grains (siltstones and shales) and carbonate lithics (largely micritic limestones and calcite spars). Shales are well lithified, well rounded, with an evident iso-orientation of the clay minerals, and for these characters, they appear to have a detrital origin, derived from the erosion of the pelitic successions of various age cropping out in the northern Apennines. Serpentinite and volcanite grains are minor components. These petrofacies characterize the western area, dominated by the Secchia River paleochannels, near Mirandola, Cavezzo and Concordia sulla Secchia, but indicate an input of the Panaro River in the San Felice sul Panaro area. The quartz-rich lithoarenites show a higher quartz-feldspar content and a lower content of siltstones, shales and carbonate lithics. This petrofacies typifies the samples from the eastern sector attributable to the Reno River paleochannels in the area between San Carlo and Mirabello. The quartzarenite sands are characterized by a significantly higher content of quartz and feldspar grains associated with metamorphic rock fragments, abundant micas and heavy minerals. These latter samples are attributable to the Po River paleochannels in the northern sector of the plain.

In situ tests
The great majority of the in situ tests presented in this work consist of CPTUs and CPTEs that are plotted in each sheet of the 56 analyzed sites, where in situ tests are available (supplement of this paper), in terms of: corrected cone resistance (q t ), sleeve friction (f s ), pore pressure (u 2 , only for CPTU) and soil behavior type index (I c ), according to Robertson (2004). This allowed a preliminary subsoil interpretation identifying the depth of the in situ ground water table, defined as "Depth of in situ GWT" in Table 1, and the "sand-like" layers associated with I c ≤ 2.6 (Idriss and Boulanger 2008).
As an example, Fig. 7 provides the CPTU results collected at the WP63 test site, located in San Carlo (municipality of Terre del Reno). The upper 6 m depth are composed by silty sands and sandy silts (approximately I c ≤ 2.6), while a thick silty-clayey layer is encountered up to 22 m depth (I c > 2.6). A sandy body is then detected between 23 and 27 m, and finally clayey deposits can be found up to 30 m. This information, coupled with the u 2 profile, makes it possible to interpret the depth of the in situ GWT, which is approximately one meter below the ground surface. Moreover, when available the sheets for each area (see supplement of this paper) include the DMTs and SDMTs geotechnical parameters according to Marchetti et al. (2001), including: the material index (I D ), the constrained modulus (M), the undrained shear strength (s u ), the horizontal stress index (K D ) and the shear wave velocity (V S ) from seismic dilatometer (SDMT) testing. Finally, when boreholes and standard penetration tests (SPT) are available, each sheet also plots an interpreted borehole log and the SPT blow count (N SPT ). Altogether, this information supports the geotechnical characterization of the liquefied sites. The compositional fields of modern and paleoriver of the Emilia plain are also reported (Lugli et al. 2007). Q: quartz; F: feldspars; L: siliciclastic rock fragments; C: carbonate rock fragments

Liquefaction assessment
At each liquefied test site, where in situ tests were available and lateral spreading was not observed (45 sites, see Table 1), liquefaction susceptibility analyses were carried out using CPT-or DMT-based methods derived from the Seed and Idriss (1971) simplified procedure with the final goal of detecting the liquefied layer(s) for the 2012 earthquakes. The test sites located approximately east of San Felice sul Panaro were associated with the mainshock that occurred on the 20 th May 2012, while the liquefied sites west of San Felice sul Panaro were coupled with the 29 th May 2012 earthquake (Pizzi and Scisciani 2012).
For the definition of the peak ground acceleration (a max ) we used the most updated ground motion model for shallow crustal earthquakes in Italy (ITA18, Lanzano et al. 2019). This ground-motion prediction equation accounts for the more recent earthquakes that have occurred in Italy, including the 2012 Emilia earthquakes, and it is calibrated by regression of empirical ground-motion amplitudes against a set of predictor variables such as earthquake magnitude, focal mechanism, source-to-site distances, and local soil conditions (i.e. the time-averaged shear wave velocity to 30 m depth). Figure 8 compares the a max from ITA18 with both the recorded data, acquired during the May 20 th and the May 29 th mainshocks, and the computed ShakeMaps (http:// shake map. rm. ingv. it/ shake4, last access Oct. 22th, 2021; Michelini et al. 2020). The lack of seismic stations in the 2012 epicentral area during the 20 th May 2012 mainshock (only one permanent station was located in the area, in the municipality of Mirandola; Cultrera et al. 2014) resulted in unreliable ShakeMap estimates in the near source region (Fig. 8a), as already noted by Cultera et al. (2014). The 29 th May ShakeMap values are better constrained because of the large number of temporary seismic stations installed in the epicentral area (Fig. 8b), and they range within one standard deviation of the ground motion model prediction. We are then confident about using the ITA18 prediction for both earthquakes, with a careful choice of the input parameters, such as the M w and the distance metrics.
Moment magnitude (M w ) of 6.1 and 5.9 are computed by Pondrelli et al. (2012) by means of the Regional Centroid-Moment Tensors (RCMT, http:// rcmt2. bo. ingv. it/, last access: October 2021; Pondrelli 2002). The M w estimates from RCMT are equivalent to  Ekström et al. 2012), which is considered the most authoritative agency for M w of earthquakes worldwide and whose magnitudes are used as a reference in many seismological studies (Di Giacomo et al. 2021). Moreover, both estimates almost overlap for M w > 5.4 (Gasperini et al. 2012).
Regarding the distance metrics, we used the Joyner-Boore distance (Rjb) to account for the effect of the extended source in ITA18, as the sites analyzed in this study are very close to the seismic sources of the mainshocks. The fault geometries, reported in Fig. 1, are taken from the available literature and limited to the area with source slip greater than 0 m: Pezzo et al. (2018) for the M w 6.1 on May 20 th and Paolucci et al. (2015) for the M w 5.9 on May 29 th .
As a consequence, the cyclic stress ratio for a M w 7.5 earthquake (CSR 7.5 ) was evaluated using the previously noted mainshock values of M w and a max . The present research is limited to reproducing the liquefaction evidence induced by the mainshoks of May 20 th and May 29 th using the simplified procedure. Reproduction of the possible additional liquefaction effects that occurred in conjunction with the aftershocks within the first four minutes after the mainshock of May 20 th would have necessitated numerical modelling to reproduce the entire pore pressure build up (i.e. from the starting time of the mainshock to at least few minutes after the last strong aftershock here considered). The performance of numerical modelling for the numerous sites analyzed in the manuscript would have required: (1) a large number of seismic recording stations in the epicentral area of the 20 th May 2012 mainshock; and (2) substantially more geotechnical information extending down to the bedrock. Unfortunately, these data are not available in the literature due to the lack of seismic stations (only the Mirandola station recorded the seismic event, see Cultrera et al. 2014) and due to the presence of deep bedrock (> 100 m in the anticlinal area of Mirandola and Casaglia, and >> 300 m in the synclinal areas that cover most of the epicentral area, see Martelli 2021. Additionally, the magnitude scaling factor (MSF) and the shear stress reduction coefficient (r d ) were evaluated according to the equations recommended by Idriss and Boulanger (2008) for CPT, DMT and CPT-DMT methods, while the depth of the earthquake ground  (Lanzano et al. 2019) water table, defined as "Depth of earthquake GWT" in Table 1, was assumed equal to 3 m depth in the case of a fluvial ridge and to 1 m in the case of a flat interfluvial depression, considering that the geomorphological features of the study area based on LIDAR data (Table 1).
Finally, the cyclic resistance ratio (CRR) for a M w 7.5 earthquake (CRR 7.5 ) was evaluated using (1) the normalized overburden corrected cone tip resistance q c1N calculated from the Idriss and Boulanger (2008) CPT-based method; (2) the horizontal stress index K D estimated from Marchetti (2016) with DMT methods; and (3) the combination of q c1N and K D parameters in the Marchetti (2016) CPT-DMT correlation. The ratio between CRR 7.5 and CSR 7.5 defines the safety factor against liquefaction (FS liq ). To screen out "clay-like" soils, a threshold was set at I c ≤ 2.6 for CPT data and at I D ≥ 1.0 for DMT measurements. Liquefaction severity indexes were also calculated for CPT profiles in terms of: • the liquefaction potential index (LPI) proposed by Iwasaki et al. (1982), given by the equation: , and z is the depth below the ground surface; • the liquefaction induced vertical settlements (S) according to Zhang et al. (2002), given by the equation: where ε vi is the post liquefaction volumetric strain, calculated as a function of the liquefaction safety factor (FS liq ) and of the equivalent clean sand normalized CPT penetration resistance (q c1N ) cs , ∆z i is the thickness of the sublayer i, and j is the number of soil sublayers; • the liquefaction severity number (LSN) proposed by van Ballegooy et al. (2014), given by the equation: where ε v is the post liquefaction volumetric strain calculated according to Zhang et al. (2002). The predictive indexes were computed to verify how well they fit the experimental evidence of the 2012 Emilia seismic sequence. As an example, Fig. 9 reports the results of the CPTU liquefaction analysis performed at the WP63 San Carlo site for the seismic event of the 20 th May 2012 (M w = 6.1). A surface non-liquefiable layer with a thickness H 1 = 1.0 m and a liquefiable layer of a thickness H 2 = 5.3 m (depth range: 1.0-6.3 m) are detected with associated LPI = 23.56, LSN = 42.39 and S = 0.27 m. By looking at the observed liquefaction evidences through the identified liquefaction-induced land damage categories (see Table 1), the liquefaction severity indexes (LPI, LSN and S) results in reasonable agreement with the 2012 effects on ground surface, classified in Category 4 (i.e. no lateral spreading but large quantities of ejected material). Reasonable agreement is noted for approximately the 75% of the analyzed cases, in terms of LPI and S, while for the remaining percentage the liquefaction prediction results underestimated or overestimated. The LSN provides reasonable agreement with the observed performance for less than the 40% of the analyzed liquefied sites, while for the majority of cases this index underpredicted the field observations. This is likely because the LSN threshold values corresponding to different severities of liquefaction are based on the New Zealand data. Further details can be detected in Table 1 where the text formatting is used to describe the accuracy of the LPI, LSN and S approaches relative to the observed damage, namely: bold-italics is underpredicted, italics is overpredicted, bold is reasonable.

Liquefaction features versus depositional facies and sediment provenance
The composition of ejected sand is a fundamental tool for the source layer identification in a sedimentary sequence, and this is pivotal for the recognition of potential areas prone to hazardous sand liquefaction phenomena. On the other hand, the petrographic composition does not necessarily appear to be a significant factor influencing the liquefaction susceptibility, as demonstrated by the wide compositional spectrum of the Po plain alluvial systems illustrated in this study. As noted previously, the stratigraphic architecture of the study area is strongly influenced by the interaction between the Po River and the network of transversal tributaries draining the Apennine chain (Fig. 10). The post-glacial Po River deposits consist of elongated channel sand bodies, with quartzarenitic composition representing detritus from wide sectors of the Alpine and Apennine chains, transported over a long distance. The Apennine sands are confined in narrow channel body, and are finer grained, containing a larger amount of silt. They show lower quartz-feldspar contents and abundant sedimentary lithics derived from relatively smaller catchments confined to the external belts of the chain; detritus is derived only from the sedimentary cover. A distinctive sand composition characterizes each river mainly based on the prevailing lithic grains. Sediments eroded from shale terrains with carbonate blocks (the chaotic Ligurian units) dominate the Secchia and Panaro catchments (litharenitic petrofacies). These lithic components are less abundant in the Reno River (quartz-rich lithoarenite petrofacies). Sand petrography data demonstrate that fluvial sand compositions have slightly varied during the Holocene with respect to the present-day sands, as shown by the compositional fields of Fig. 6. This is particularly evident for the modern Reno River sands, which are more litharenitic compared to its paleochannel. The Po River sands are richer in quartz, feldspar, and metamorphic grains.
The petrography of the sand blows therefore makes it possible to discriminate the sedimentary provenance of the sand ejected during the 2012 seismic crisis and to ascribe them to their respective Po, Secchia and Reno Holocene paleochannels (Fig. 10). The petrography of sands has also supported the identification of the source level of liquefaction, in all sites located at depth between 5 and 8 m in the alluvial sediment sequence.
This comprehensive dataset shows that earthquake-induced liquefaction phenomena affected sand layers with significant contents of non-plastic silt and was not limited to clean sands and well-sorted deposits, as suggested by previous results from Fontana et al. (2015Fontana et al. ( , 2019 and Amoroso et al. (2020). Within the considered alluvial plain, liquefaction primarily affects silty sand with litharenitic composition from laterally confined depositional bodies, such as levee, crevasse splay, to minor channel ribbons. Subordinately liquefaction was related to thick channel-fill sands deposited by the Po river.  Stefani et al. 2018). Liquefaction sites are aligned along paleochannels as also reported by Civico et al.(2015) 6 Liquefaction charts The liquefied layers identified at the trial sites (Table 1) were subsequently used to develop summary liquefaction charts, as reported in Figs. 11 and 12. In particular, Fig. 11a compares the average CSR 7.5 (CSR 7.5 avg ) and average q c1N (q c1N avg ) within each identified liquefied layer, with the CPT-based liquefaction triggering curves by Idriss and Boulanger (2008). The average parameters (i.e. CSR 7.5 avg and q c1N avg ) were calculated using  Idriss and Boulanger (2008), (b) fines content correlations by Suzuki et al. (1998) and Boulanger and Idriss (2014). Dashed lines indicates the upper and lower boundaries of the two curves  (2008) triggering curves which indicate liquefaction in all cases, even though the plot provides CSR 7.5 avg and q c1N avg pairs for the entire liquefied layer at each site, rather than the "critical layer" computed with the lowest average q c1Ncs /CSR 7.5 ratio (where q c1Ncs is the normalized cone tip resistance for clean sand) over a depth range of about a meter. The identification of the critical layer at each site would have produced q c1N avg values somewhat smaller (less than 38%) with little effect on the CSR 7.5 avg (less than 9%), but the results would still plot above the liquefaction triggering curves for all the data points indicating agreement with the observed behavior. Figure 11b plots the FC from the sand boil samples versus the average I c (I c avg ) estimated from the 2012 liquefied layer in comparison with correlations by Suzuki et al. (1998) and by Boulanger and Idriss (2014). The FC and I c avg data pairs exhibit considerable scatter for a significant number of case studies relative to the average curves provided by Suzuki et al. (1998) and Boulanger and Idriss (2014) (Fig. 11b), although less than the 35-40% of the cases plot outside the upper and lower boundary curves. This scatter may be related to a site-specific dependency of these kind of I c -FC correlations (Idriss and Boulanger 2008;Boulanger and Idriss 2014) and to the fact that I c can be seen as a parameter related to the mechanical soil response, including soil plasticity, and not strictly to the grain size of the soil deposits alone (Robertson and Cabal 2015).
Finally, Fig. 12 evaluates liquefaction manifestation at the ground surface relative to the presence of a non-liquefiable surface crust if thickness H 1 overlying a liquefiable deposit of thickness H 2 according to charts developed by Ishihara (1985) and by Towhata et al. (2016) using H 1 and LPI. Almost the 90% of the experimental liquefaction observations correctly fit with Ishihara chart (Fig. 12a), considering the range of peak ground accelerations used in the liquefaction assessment (a max = 0.40-0.46g, see also Table 1). In only four cases, related to the municipalities of Bondeno, Cavezzo and Terre del Reno, does the chart incorrectly predict that liquefaction would not be observed at the ground surface probably due to the presence of a thick non-liquefiable crust (H 1 ≈ 4-10 m). The prediction provided by the Towhata chart (Fig. 12b) is far less accurate and contemplates highly probable surface manifestations of liquefaction for less than 25% of the analyzed cases. In this respect, for the analyzed Emilia dataset the chart proposed by Ishihara (1985) looks to provide more reliable results using the H 1 -H 2 pairs than the one proposed by Towhata et al. (2016) where the non-liquefiable crust is coupled with the LPI in place of H 2 .

Conclusions
We have compiled a large and comprehensive dataset consisting of 120 liquefied sites induced by the 2012 Emilia Romagna earthquakes in the alluvial plain created by the late Quaternary sedimentation of Apennine fluvial systems and of the Po River. This compilation includes the entire available structural, grain size, compositional and geotechnical data for each liquefaction site consisting of a wide spectrum of sediment types.
A strong relationship exists between the occurrence of liquefaction and the stratigraphic architecture of the subsurface. The granular sediments ejected to the surface are entirely derived from Holocene units, associated with river channel, levee, and crevasse deposits. The geometry of the fluvial paleochannels contributed to confine the liquefaction sites along narrow ribbons crossing the plain. In addition, coseismic surface fractures were largely localized across fluvial ridge morphologies. The source level for liquefaction in all sites was identified at shallow depth, in sand layers, mainly between 3 and 6 m, covered by fine-grained cohesive deposits. The granular sediments involved in the liquefaction phenomena often show a significant amount of non-plastic silt. The texture and composition of sand blows has allowed us to identify the provenance of the sand ejected and to ascribe them to their respective fluvial channel. The sands deposited by the Po River are organized in large, vertically stacked channel bodies with quartzarenitic composition, whereas the Apennine drainage systems generated narrower channel bodies rich in finer grained litharenitic sands and silty sands.
The availability of in situ tests at the liquefied sites allows us to verify and validate the reliability of the liquefaction charts in alluvial sediments. In this respect, the comparisons of the liquefaction severity indexes with the coseismic observations revealed that the LPI and S generally is consistent with the 2012 evidences while the LSN seems to underestimate the surface manifestation of liquefaction probably due to the dependency of LSN on the New Zealand dataset for which it has been developed.
The CPT data points at the liquefied sites fit well with the triggering chart developed by Idriss and Boulanger (2008), highlighting the consistency of this simplified method for evaluating liquefaction. In contrast, the I c -FC correlations by both Suzuki et al. (1998) and Boulanger and Idriss (2014) only provide a broad estimate of the true fines contents of the 2012 sand ejecta (as measured by the grain size analyses), although 60-65% of measured fines content from sand boils at liquefaction sites fall within the upper and lower boundary curves. The dispersion of these data points around the average curves denotes the site-specific nature of I c -FC correlations and the effect of plasticity and mechanical behavior associated with I c rather than simply the grain size.
For the vast majority of the data points (almost the 90% of the analyzed cases) for the Emilia earthquakes, the presence of a medium-thick non-liquefiable crust (≈ 3-6 m thick) related to the range of the 2012 peak ground accelerations (a max = 0.40-0.46g) did not prevent surface manifestations of liquefaction in agreement with the liquefaction charts proposed by Ishihara (1985). In contrast, the Towhata et al. (2016) chart only predicted surface liquefaction features for a small percentage of the dataset (less than 25% of the analyzed cases) where liquefaction actually occurred. Therefore, for the analyzed Emilia case studies the use of non-liquefiable crust provides better estimations of the liquefaction manifestations when coupled with the thickness of the liquefiable layer (H 2 ) rather than with the liquefaction potential index (LPI). These coseismic liquefaction observations emphasize the importance of an in-depth study of geological and geotechnical properties of these crusts.