Liquefaction susceptibility assessment in fluvial plains using airborne lidar: the case of the 2012 Emilia earthquake sequence area (Italy)

We report a case study from the Po River plain region (northern Italy), where significant liquefaction-related land and property damage occurred during the 2012 Emilia seismic sequence. We took advantage of a 1 m pixel lidar digital terrain model (DTM) and of the 2012 Emilia coseismic liquefaction data set to (a) perform a detailed geomorphological study of the Po River plain area and (b) quantitatively define the liquefaction susceptibility of the geomorphologic features that experienced different abundance of liquefaction. One main finding is that linear topographic highs of fluvial origin – together with crevasse splays, abandoned riverbeds and very young land reclamation areas – acted as a preferential location for the occurrence of liquefaction phenomena. Moreover, we quantitatively defined a hierarchy in terms of liquefaction susceptibility for an ideal fluvial environment. We observed that a very high liquefaction susceptibility is found in coincidence with fluvial landforms, a high-tomoderate liquefaction susceptibility within a buffer distance of 100 and 200 m from mapped fluvial landforms and a low liquefaction susceptibility outside fluvial landforms and relative buffer areas. Lidar data allowed a significant improvement in mapping with respect to conventionally available topographic data and/or aerial imagery. These results have significant implications for accurate hazard and risk assessment as well as for land-use planning. We propose a simple geomorphological approach for liquefaction susceptibility estimation. Our findings can be applied to areas beyond Emilia that are characterized by similar fluvial-dominated environments and prone to significant seismic hazard.


Introduction
Soil liquefaction is one of the most outstanding hydrogeological processes that can occur during earthquakes, provided the existence of saturated loose sandy layers confined by impermeable deposits within a certain distance from the earthquake epicenter.In fact, it is well known that seismic shaking during earthquakes can cause water-saturated sediments (low density or uncompacted sandy and/or silty sandy deposits) to temporarily lose strength and to act as a fluid.Because of the peculiar subsurface stratigraphy, this is a common occurrence in alluvial and coastal plains (De Martini et al., 2012).
The assessment of liquefaction susceptibility in fluvial and coastal plains is a worldwide topic that is growing in interest, probably because of the recent need for "new" areas for urban and industrial development related to the fast growing of global population.The scientific literature clearly highlights the great importance of studies on liquefaction as a contribution to seismic hazard assessment at the local and regional scale (for instance the epicentral area of historical or modern earthquakes - Kotoda et al., 1988;Obermeier, 1996).In fact, geotechnical/stratigraphical and geomorphological studies in liquefaction-prone areas can provide critical information that will be potentially useful for making accurate hazard and risk maps for land-use planning.
Several approaches have been proposed so far in order to predict the liquefaction susceptibility at a specific site and to elaborate liquefaction hazard maps.From a geotechnical point of view, the preferred ones are the liquefaction potential index (LPI) and the liquefaction severity number (LSN) (Papathanassiou et al., 2015, and references therein).Both indexes are based on the geotechnical engineering properties of soil derived from methods such as the standard penetration test (SPT) and the cone penetration test (CPT).Differently, few studies have applied a mainly geomorphological approach for the identification of liquefaction-prone regions.Wakamatsu (1992) grouped sedimentary deposits into three categories of liquefaction susceptibility (likely, possible and not likely) using geomorphological criteria and classifying geomorphological units such as natural levee, former river channel, sandy dry river channel and artificial fill as the highest level of liquefaction potential.Witter et al. (2006) and Ganapathy and Rajawat (2012) combined geomorphic expression and geological characteristics to define and distinguish lithological units and to produce liquefaction susceptibility maps for the central San Francisco Bay region (USA) and Chennai city (India), respectively.
Historical accounts, maps and aerial photographs allowed Wotherspoon et al. (2012) to identify areas of land reclamation and old channels that have had flow diverted away in and around the town of Kaiapoi (north of Christchurch, New Zealand).Their findings show that these areas correlated well with many of the areas having significant liquefaction damage following the 2010 Darfield earthquake (New Zealand).
In this paper we report a case study from the Po River plain region (northern Italy), where significant liquefactionrelated land and property damage occurred during the 2012 Emilia seismic sequence.This case highlighted the need for a better understanding of the liquefaction hazard, especially in complex fluvial environments (EMERGEO Working Group, 2013, and references therein).Notably, the study area is dominated by agricultural fields, and thus the majority of the geotechnical data are present only in small towns and villages, preventing a comprehensive estimate of liquefaction potential over the entire 2012 coseismic area.
The huge amount of data on coseismic liquefaction related to the May-June 2012 events (EMERGEO Working Group, 2012a, 2013) offers a unique opportunity to refine our knowledge and methodologies to better understand how the geomorphological setting, besides the stratigraphic and hy-drologic conditions, can favor the occurrence of liquefaction phenomena.
Thus, using a 1 m pixel airborne lidar data set, we analyzed in detail the correlation between geomorphological features and liquefaction occurrence in the area affected by the 2012 Emilia seismic sequence with the following objectives: (a) perform a detailed geomorphological study of the Po River plain sector that experienced liquefaction phenomena; (b) quantitatively define the areas where liquefaction is more likely to occur in terms of liquefaction susceptibility hierarchy.
In addition, we introduced a simple geomorphological approach for the estimation of liquefaction susceptibility using the liquefaction density (Ld -observed liquefaction effects over a certain area) since it can be calculated entirely by means of remote sensing techniques, to be combined with in situ LPI/LNS indexes calculation, where available.
This paper is organized as follows: Sect. 2 describes the geomorphological setting of the epicentral area and provides the background to the 2012 Emilia earthquake sequence and the details of the observed liquefaction phenomena; Sect. 3 presents the geomorphological analysis, including landform mapping based on a high-resolution lidar digital terrain model (DTM); and Sect. 4 presents the results of the quantitative analysis of the spatial distribution of the 2012 coseismic observations.The paper closes with a discussion of the main findings of this work.

Geomorphological setting
The area struck by the 2012 seismic sequence is characterized by an alluvial plain with flat morphology (average slope less than 3 • ) and by drainage and fluvial landforms related to the Po, Secchia, Panaro and Reno rivers (Regione Emilia-Romagna, 1999).This landscape has been strongly influenced by human activity since ancient Roman times, by the drainage and reclamation of swamps and marshes, as well as by several diversions of river courses and by the creation of artificial canals and levees.
In more detail, the study area is dominated by a complex drainage and paleo-drainage pattern, characterized by the presence of long and smooth alluvial ridges rising over the floodplain, abandoned riverbeds, levees (both natural and artificial) and crevasse splays.Alluvial ridges are the result of the action of aggrading rivers and are characterized by the rising of the streambed elevation due to deposition of sediments, whereas abandoned riverbeds are the result of a process known as avulsion, which is a shift in the main channel bed of the river to a new course.Levees are wedge-shaped ridges of sediment bordering river channels, and they generally occur as sinuous, ribbon-like prismatic bodies.They are among the coarser sediments because of their position proximal to the active channel.Crevasse splays form when the levee of the channel is breached and water carrying sediments flows out onto the floodplain, depositing sediments in a similar pattern to an alluvial fan deposit.
The general geological-stratigraphic setting of the area is characterized by a subsoil consisting of Holocene and poorly consolidated fluvial deposits, getting generally younger from west to east.In more detail, the uppermost tens of meters of subsoil is composed by clay, silt and sand deposits organized in layers and lenses, related to the frequent migration of the rivers in the alluvial plain (Regione Emilia-Romagna, 1999;De Martini et al., 2014).Predominantly sandy and silty-sand deposits are generally found in coincidence with river channels, levees and crevasse splays, whereas all the areas dominated by clayey or silty-clayey deposits are related to lowenergy environments quite often "trapped" between active or paleo-riverbeds.
In general, it is possible to identify two overlying aquifer systems in the study area.The shallower one is a semiconfined and locally phreatic aquifer consisting of interconnected silty-sand lenses with a variable thickness up to 6-8 m.It is underlain by an aquiclude composed of silts, clays, and peats, which confine the second, deeper aquifer, located at a mean depth of 16-18 m b.g.l.(Papathanassiou et al., 2012).

The 2012 Emilia seismic sequence and observed liquefaction phenomena
In May-June 2012, a seismic sequence struck a broad area of the Emilia-Romagna region in northern Italy (Fig. 1), resulting in 26 fatalities and hundreds of injured, 15 000 homeless, severe damage of historical centers and industrial areas, and an estimated economic toll of more than EUR 2 billion.The sequence included two mainshocks ML 5.9 and ML 5.8, on 20 and 29 May, respectively) and five ML ≥ 5.1 events.The maximum acceleration for the 20 May earthquake was recorded in the town of Mirandola (located ca. 13 km from the epicenter), with peak ground acceleration of 0.31 and 0.26 g along the vertical and horizontal components, respectively (Bozzoni et al., 2012).The whole aftershock area, as defined by more than 1800 ML > 1.5 earthquakes (ML after Mazza et al., 2012), extended over more than 50 km, being elongated in the WNW-ESE direction.The focal mechanisms of the main earthquakes consistently show a compressional kinematics with E-W-oriented nodal planes (Pondrelli et al., 2012;Scognamiglio et al., 2012;TDMT Database, 2015).These mechanisms are consistent with a horizontal, N-S-oriented compression (P axes) defined also by presentday stress indicators in the region (Montone et al., 2012, and references therein) and by the GPS-derived velocity field (Devoti et al., 2011).The 2012 seismic sequence confirmed the Po Plain to be prone to coseismic liquefaction phenomena, as also mentioned in several historical reports of earthquakes that have occurred in northern Italy, such as the Ferrara 1570, Soncino 1802 and Salò 1901 events (Galli, 2000).
The 2012 Emilia earthquake sequence triggered widespread liquefaction through a vast area of the Po River plain (Emilia-Romagna and Lombardia regions).An area of approximately 1200 km 2 was covered by the 2012 post-earthquake surveys, and coseismic geological evidence was collected through field reconnaissance and aerial surveys, reports from local people and Web-based surveys (EMERGEO Working Group, 2012a, 2013).More than 1350 coseismic geological effects were identified at the surface (Fig. 1).These consist of liquefaction phenomena induced by seismic shaking.
On the basis of their morphologic and structural characteristics, the observed coseismic effects at the surface (Fig. 2) can be grouped into three main classes (for details see The surveyed features appear independent from the type of environment, as they occur on roads, buildings, backyards, parks, agricultural fields, etc.Some manmade underground structures -such as wells, foundations, and sewers -forming artificial boundaries represent a simpler escape for the overpressured water and sediments. Thanks to the information provided by several local eyewitnesses, collected by the EMERGEO personnel in the field, and based on the magnitude and timing of the main earth-quakes, we know that the liquefaction process was induced by the 20 and 29 May mainshocks only.In fact, the only large aftershock that could have produced liquefaction is the 20 May 2012, ML 5.1, event that occurred only 4 min after the first mainshock, being clearly indistinguishable to eyewitnesses.In addition, we have evidence for only five sites where sand blows reactivated following the second mainshock (29 May).Liquefaction reactivations occurred only in the San Felice sul Panaro municipality (see Fig. 4 in EMERGEO Working Group, 2013), which is located less than 10 km from both mainshocks.
Large quantities of ejected sand, silt and water produced damage to commercial buildings, residential houses and infrastructure within the urbanized area.In more detail, subsidence, uplift and lateral spread due to liquefaction of soil at relatively shallow depth were found.Most of the observed damage involved masonry buildings, precast structures, and in some cases reinforced concrete buildings.The overall damage to residential buildings was minor, but damage was particularly evident for old and poorly maintained masonry structures.On the other hand, significant damage to industrial facilities was observed (Cimellaro et al., 2014).
It was already qualitatively noticed that the presence of subtle fluvial landforms in the area hit by the 2012 seismic sequence strongly influenced and favored the occurrence of liquefactions (Bertolini and Fioroni, 2012;Di Manna et al., 2012;Ninfo et al., 2012;De Martini et al., 2014;EMER-GEO Working Group, 2013;Papathanassiou et al., 2012).In fact, looking at the distribution of the liquefaction features, it appears clearly that they are not evenly distributed over the area but are mostly arranged in clusters and rectilinear or meander-like alignments (Fig. 3).
Moreover, the maximum distance of the observed liquefaction from the earthquake(s) epicenter is about 30 km.This is in good agreement with national and international empirical relationships available (Obermeier, 1996;Galli, 2000) between earthquake magnitude and maximum distance for the occurrence of liquefaction phenomena.

Methodology
The non-homogeneous spatial distribution of the 2012 liquefaction phenomena (see Sect. 2.2) in an area of similar stratigraphy and hydrologic conditions highlighted the need for a better understanding of the factors favoring liquefactions.Therefore, we developed a detailed geomorphological characterization of the area, with special focus on subtle landforms of fluvial origin.
We performed a substantial refinement in mapping of the fluvial landforms already identified by Castiglioni et al. (1999), by using a high-resolution (1 m) airborne lidar topography, and defined the geomorphologic characteristics of those areas that experienced liquefaction phenomena during the 2012 Emilia seismic sequence.

High-resolution topography analysis
Widely available digital elevation data (e.g., 10 m digital elevation model -DEM of Italy; Tarquini et al., 2007) are generally too coarse to allow identification and precise mapping of fluvial landforms in alluvial plains (characterized by low morphological gradients).To overcome this drawback, we took advantage of a 1 m resolution lidar DTM.The airborne lidar survey of the study area was performed in 2008 by Regione Emilia-Romagna using an Optech Gemini ALTM sensor and covered an area of 693 km 2 (lidar accuracy: 0.15 m vertical, 0.35 m horizontal).We processed the original 1 × 1 m pixel DTM in order to obtain several derivative digital maps (shaded relief, slope, aspect, etc.) useful to reveal even the subtlest (0.5 m in elevation) topographic features.
In addition, we computed a residual DTM (residual relief surface -RRS -Fig.4) by removing the regional orographic trend.We eliminated the large-scale variations from the original DTM by calculating a smoothed surface (Grohman and Riccomini, 2009), and we then derived a residual DTM subtracting the smoothed DTM from the original DTM.The residual DTM obtained (RRS) emphasizes the fluvial landforms (both rises and depressions) and makes their comparison at different locations easier, independently from the regional trend.This is clearly visible in Fig. 4, with the RRS showing in a better way depressions and rises over the floodplain and helping substantially the visual interpretation.adopted by the previously published 1 : 250 000 geomorphological map of the Po Plain (Castiglioni et al., 1999), but it substantially integrates and, in some cases, deeply revises it.

Identification and mapping of fluvial landforms
In detail, according to Castiglioni et al. (1999), landforms in the study area are classified as fluvial and fluvioglacial, aeolian, tectonic and anthropic.We focused our efforts on the precise identification and detailed mapping of selected fluvial landforms either rising over the surrounding floodplain (levee ridges and alluvial ridges, crevasse splays) or incised (traces of abandoned riverbeds).
The map we present comprises an area of 693 km 2 that was surveyed at 1 : 25 000 scale, as the best compromise between the available lidar data set, the desired detail and the time needed to survey the investigated area.We performed mapping mainly based on lidar-derived topography and aerialphoto interpretation.We adopted a two-step approach: we first utilized the residual topography to obtain a first-order identification of concavities (depressed landforms) and convexities (rising landforms) in the area, and then we precisely mapped on screen the selected fluvial landforms.Figure 5 is an example of the improved identification and mapping of landforms thanks to high-resolution (1 m) lidarderived DTM.Subtle (i.e., less than 1 m relief) features like abandoned riverbeds (Fig. 5b-c) are clearly revealed in the shaded relief image (Fig. 5a-b) as well as by the topographic profile (Fig. 5c).
We mapped in detail more than 100 km 2 of fluvial landforms.The result of our analysis is shown in the map of Fig. 6, together with the observed 2012 liquefaction phenomena (EMERGEO Working Group, 2013).

Results
In order to find a simple and homogeneous approach to define areas characterized by different liquefaction potential, we performed a quantitative analysis of the spatial distribution of the 2012 coseismic liquefaction effects using GIS tools.We analyzed the spatial relationships (overlay and proximity analysis) between our precisely located fluvial landforms and the location of the 2012 liquefaction effects (Fig. 6), identified from both field and aerial surveys, as published by the EMERGEO Working Group (2013).Taking into account the areal extent of the available lidar data, we made a subset from the EMERGEO Working Group ( 2013) data set of 1306 coseismic effects out of a total of 1350.
As a first consideration, the analysis of the spatial distribution of the liquefaction effects shows that 699 out of a total of 1306 observed liquefaction phenomena (53 %) are located exactly in coincidence with mapped fluvial landforms, which notably represent only ∼ 15 % of the whole study area.
Among the liquefactions observed on mapped fluvial landforms, alluvial ridges and levee ridges hosted 63 % of the observed liquefaction effects, while crevasse splays account for 20 % and abandoned riverbeds for 17 % (Fig. 7).
As for the liquefaction effects observed outside mapped fluvial landforms, most of them (about 500) appear randomly distributed over the floodplain.Conversely, less than 100 liquefaction effects show a spatial distribution (e.g., meander-like alignments) that can potentially be related to concealed/undiscovered fluvial features.
Furthermore, in order to homogeneously describe the relationship between 2012 liquefaction phenomena and fluvial landforms and to find a simple approach for the assessment of liquefaction susceptibility to be adopted for future earthquakes, we decided to compute a Ld parameter, set as the number of observed coseismic liquefaction effects over a certain area.In general, the whole investigated area (∼ 693 km 2 ) shows an average Ld of 1.9 effects km −2 (1306 effects 693.82 km −2 ).We then computed the Ld of the mapped fluvial landforms, obtaining a value of 6.6 ef- fects km −2 (699 effects 106.04 km −2 ), 3.5 times larger than the value computed for the whole area.
Comparing these two values, it is possible to affirm that specific fluvial landforms -such as abandoned riverbeds, levee ridges and alluvial ridges and crevasse splays -are a preferential location for the occurrence of liquefaction vents.
In addition, we set two arbitrary buffer distances from the fluvial landforms (100 and 200 m) in order to investigate the behavior of the areas located in close proximity to mapped fluvial landforms, and thus we computed the relative Ld.A Ld of 2.3 effects km −2 characterizes the buffer distance of 100m (133 effects 57.62 km −2 ), whereas a Ld of 1.7 effects km −2 is obtained within a distance of 200m from fluvial landforms (96 effects 55.96 km −2 ), these values being 2.9 and 3.9 times smaller than the Ld computed for fluvial features, respectively.Finally, we noticed that, when computed outside the 200 m buffer area, the Ld drops down to a value of 0.8 effects km −2 , notably 1 order of magnitude smaller than the value obtained in coincidence with specific fluvial landforms.

Discussion
In this paper, we investigated the correlation between geomorphological features and liquefaction occurrence in the Po River plain, which was affected by the 2012 Emilia seismic sequence.In doing this, we took advantage of the huge amount of data on coseismic liquefaction effects related to the 2012 Emilia seismic sequence.Moreover, we put special emphasis on the fluvial features, already recognized as a potential control factor for the preferential location of liquefaction phenomena (Bertolini and Fioroni, 2012;Di Manna et al., 2012;Papathanassiou et al., 2012;EMERGEO Working Group, 2013;De Martini et al., 2014).Ridges, crevasse splays and riverbeds are usually richer in water-saturated coarse layers with respect to the whole flat alluvial plain (dominated by overflooding unsaturated fine deposits), and for this reason these fluvial features are commonly prone to liquefaction.
We quantitatively defined the liquefaction susceptibility of the geomorphologic features that experienced different abundance of liquefaction using the Ld.In fact, a hierarchy in terms of liquefaction susceptibility could be set for an ideal fluvial environment by using the ratio between any Ld and the maximum Ld as calculated for the 2012 Emilia case study.We set three classes of susceptibility: very high, high to moderate, and low, where the ratio is 1 for very high, between 0.99 and 0.15 for high to moderate, and ≤ 0.15 for low. Figure 8 shows an example of this liquefaction susceptibility assessment from the Sant'Agostino area (southeast of the study area), where liquefaction phenomena were particularly abundant and severe, this channel being a very young land reclamation area (beginning of the nineteenth century).Here, Lds are 19.4 for specific fluvial landforms (alluvial and levee ridges, crevasse splays), 16.7 within a buffer distance of 100 m, 8 within a buffer distance of 200 m and 2 outside mapped fluvial landforms and relative buffers.These densities translate to very high susceptibility (ratio 1) for fluvial landforms, high to moderate susceptibility (ratio 0.8 to 0.4) for buffer 100 and 200 m, and low susceptibility (ratio 0.1) outside fluvial features and buffers.
The relation between old and/or former river channels/landforms and liquefaction occurrence was already recognized for some historical and modern earthquakes.Liquefaction-induced damage was observed, among others, during the 1990 Luzon, Philippines earthquake (M7.8); the 2007 Niigata-Chuetsu-Oki, Japan earthquake (M6.8;Orense et al., 1991Orense et al., , 2008)); and the 2010-2011 New Zealand earthquakes (Wotherspoon et al., 2012;Villamor et al., 2014), showing that coastal and alluvial plains are morphological settings favorable to earthquake-induced ground cracks and liquefaction phenomena.Furthermore, coastal and allu-vial plains are often a preferred location for industrial/power plants and urban centers, thus highlighting the importance of knowing the location of old and/or former channels and fluvial landforms when defining liquefaction-prone areas.
We present a simple geomorphological approach for the estimate of liquefaction susceptibility.DTMs derived from lidar are going to greatly improve our ability to better identify and map landforms that are prone to liquefaction, and we may need to use high-resolution geomorphic mapping either as a first/preliminary approach to estimate liquefaction susceptibility (especially in areas with poor or no geotechnical data) or as a way to refine maps that are based on geotechnically derived indexes.
Finally, we propose to test our findings in future works by comparing to already available LPI/LNS values and to export the same approach in similar fluvial-dominated environment that are prone to significant seismic hazard.

Conclusions
On the basis of the analysis of the spatial distribution of the 2012 coseismic observations, one main finding of our study is that linear topographic highs of fluvial origin (alluvial ridges and levees) -together with crevasse splays, abandoned riverbeds and very young land reclamation areas -hosted approximately half (53 %) of the total number of observed liquefaction features, which notably represent only ∼ 15 % of the whole study area, and thus they acted as a preferential location for the occurrence of liquefaction phenomena.As for the liquefaction effects observed outside mapped fluvial landforms, apart from less than 100 liquefaction effects showing a spatial distribution potentially related to undiscovered fluvial features, most of them (about 500) appear randomly distributed over the floodplain.Knowing that alluvial ridges and levees attracted the development of urban settlement and associated infrastructure since ever, and that liquefaction phenomena occurred preferentially within or in proximity to them, the related risk appears to be quite important.
Moreover, we quantitatively defined a hierarchy in terms of liquefaction susceptibility for an ideal fluvial environment by using the 2012 Emilia data.We observed that a very high liquefaction susceptibility is found in coincidence with fluvial landforms, a high to moderate liquefaction susceptibility within a buffer distance of 100 and 200 m from fluvial landforms and a low liquefaction susceptibility outside fluvial features and relative buffer areas.
Lidar data analysis substantially increased the confidence in identifying even subtle (sub-metric) fluvial geomorphic features in areas of very low relief.This resulted in a significant improvement of mapping with respect to conventionally available digital topographic data and/or aerial-photo interpretation.
We present a simple geomorphological approach for the estimate of liquefaction susceptibility.These results have sig-nificant implications that will be potentially useful for accurate hazard and risk assessment, in order to avoid or mitigate liquefaction-induced damage.Beside the Po River alluvial plain, the findings of this study can be applied to areas beyond Emilia that are characterized by similar fluvialdominated environments and prone to significant seismic hazard.
Author contributions.R. Civico carried out the lidar processing, the mapping and the spatial analysis and wrote most of the paper.C. A. Brunori contributed in lidar processing and drafted the manuscript.P. M. De Martini had the idea to analyze the correlation between geomorphological features and liquefaction occurrence and drafted the manuscript.F. R. Cinti, D. Pantosti and S. Pucci provided useful suggestions that much improved the manuscript.All the authors equally contributed to data collection in the field in the framework of the EMERGEO Working Group.
Acknowledgements.The airborne lidar survey was kindly released by Regione Emilia-Romagna.We are grateful to Fabio Villani for helpful discussions and advice.We thank NH ESS editor Oded Katz for the editorial handling, and we are grateful to Pilar Villamor for the constructive suggestions and comments that strongly improved the paper.One anonymous reviewer also provided insightful comments.The views and conclusions contained in this study are those of the authors and are purely scientific, and they should not be interpreted as necessarily representing official policies, either expressed or implied, of the Istituto Nazionale di Geofisica e Vulcanologia (INGV).The authors and the INGV make no warranty, either expressed or implied, regarding the suitability of these results for a particular use, and they shall not be liable under any circumstances for any direct, indirect, special, incidental, or consequential damages with respect to claims by users of this study.
Edited by: O.Katz Reviewed by: P. Villamor and one anonymous referee

Figure 2 .
Figure 2. Examples of observed liquefaction phenomena produced by the 2012 Emilia seismic sequence: (a, b) fractures; (c) single sand volcano; (d) alignment of sand volcanoes (photos courtesy of EMERGEO Working Group, 2012b).
Fig. 4 in EMERGEO Working Group, 2013): (a) individual sand blows, scattered vents and coalescent flat cones; (b) elongated/aligned multiple sand volcanoes, fissures with coalescing sand blows and sand flows from coseismic open fractures occurring both on natural and paved ground surfaces; and (c) newly formed open fractures and cracks without evident sand extrusion at the surface, which may be associated with subsidence, bulging or lateral spreading related to sediment liquefaction.

Figure 3 .
Figure 3. Subset of the 2012 Emilia sequence liquefaction effects (from both field and aerial surveys) showing their arrangement in clusters and rectilinear or meander-like alignments.

Figure 4 .
Figure 4. Comparison between the original 1 m lidar DTM (upper panel) and the residual relief surface (RRS -lower panel).The calculation of the RRS emphasizes rising and depressed landforms and makes easier their comparison at different locations, independently from the regional trend of the topography.

Figure 5 .
Figure 5. Example of improved high-resolution lidar-derived digital elevation model and application in geomorphological mapping from an area north of Gavello (FE), northern portion of the study area (see Fig. 1 for location): (a) 1 m lidar DTM with no interpretation; (b) same area showing subtle fluvial landforms like alluvial ridges and abandoned riverbeds; (c) topographic profile (black line in Fig. 5b) showing the subtle morphologic expression (total relief less than 1 m) of the abandoned riverbed.

Figure 6 .
Figure 6.Geomorphological map of the study area (see Fig. 1 for location) showing selected fluvial landforms and observed liquefaction phenomena.

Figure 7 .
Figure 7. Distribution of the observed coseismic liquefaction effects with respect to mapped fluvial landforms and to the floodplain (see Fig. 6).

Figure 8 .
Figure 8. Liquefaction susceptibility hierarchy map for the Sant'Agostino village area (see Fig. 1 for location).In the background the 1 m pixel DTM of the area: elevations are in gray scale from highest (white) to lowest (black) values.