A comparison between bivariate and multivariate methods to assess susceptibility to liquefaction-related coseismic surface effects in the Po Plain (Northern Italy)

ABSTRACT The two main events occurred during the 2012 seismic sequence on the Po Plain in Northern Italy (May 20 Mw 5.9 and May 29 Mw 5.8) induced widespread coseismic surface effects (CSEs), mostly liquefaction phenomena, which caused severe damages. To prevent risks related to CSE, an important non-structural measure is to assess locations susceptible to liquefaction, which is usually done at site scale by means of deterministic approaches in micro-zonation studies. This paper brings a novel methodological contribution in the field of CSE susceptibility mapping at regional scale, by testing and comparing bivariate (Weights of Evidence) and multivariate (Logistic Regression) methods that, so far, have not yet been used for such purpose. In a study area of 1480 km², the computation used an inventory of CSE as supporting evidence and a set of maps with geologic, geomorphic, hydrogeological and topographic factors as causal variables. Both methods provided susceptibility maps with a significant predictive capability and a fairly good spatial agreement between each other. In conclusion, this application of data-driven spatial modelling approaches indicates that such methods can be effectively used for liquefaction susceptibility zonation at regional scale, which can be of significant support for spatial planning over large areas.


Introduction
Coseismic surface effects (CSEs) related to liquefaction of saturated sandy deposits are among the most important indirect hazardous phenomena associated with earthquakes. In free-field conditions, CSEs include phenomena such as soil fissures, sand boils, soil bulging and lateral spreading, ejections of sand and water from wells (Huang and Yu 2013). In built environments, liquefaction causes the loss of bearing capacity of soils, thus provoking sudden settlements or even failure of buildings and severe damages to urban structures and infrastructures (Galli et al. 2012;Huang and Yu 2013;Civico et al. 2015). In general, soil liquefaction is triggered by strong earthquakes (Mw > 5.5) within a certain distance from the earthquake epicentre. In saturated sand-silt layers, located typically within the first 20 m from the ground, shaking causes rapid increase of pore-water pressure with consequent reduction of effective stress and transition from solid to liquid state (Marcuson 1978;Youd and Idriss 2001). The occurrence of CSE is therefore strongly associated to the stratigraphy of normally consolidated deposits. The presence of saturated loose silty-sand layers confined by impermeable silt-clay layers is a major predisposing condition for liquefaction.
CONTACT C. Lanfredi Sofia cinziasofia@campus.ul.pt Nevertheless, CSE occurrence in floodplain areas is also related to geomorphological features such as sandy abandoned streambeds, abandoned river levees, and buried dunes of alluvial sands (Owen and Moretti 2011;De Martini et al. 2012).
To prevent risks related to CSE, one of the most important non-structural measures is to assess and map locations susceptible to liquefaction, which is usually done, at site scale, by means of deterministic approaches for micro-zonation studies and, at regional scale, by inference of geologic and hydrogeological predisposing conditions.
The challenge of this study is to assess and map CSE susceptibility at regional scale in the Po Plain (Italy) by using bivariate and multivariate geostatistical data-driven methods. Specifically, the study area includes a 1400 km 2 portion of Emilia-Romagna Region in which, during the May 2012 seismic sequence peaking to Mw 5.9, hundreds of CSEs took place and were later on inventoried (Emergeo Working Group 2013) and also detected by satellite images (Chini et al. 2015). The availability of a georeferenced inventory of hundreds of CSEs occurred during the 2012 earthquake provided the unique opportunity to test geostatistical data-driven methods to assess liquefaction-prone areas by exploring the relationships between the location of CSE (supporting evidences) and the variability of geologic, geomorphic, hydrogeological and topographic conditions (predisposing factors). Indeed, in order to map CSE susceptibility in the study area, we have applied and compared the Weights of Evidence (WofE) (Bonham-Carter et al. 1989;Lee et al. 2002) and the Logistic Regression (LR) methods (Cox 1958;Agterberg et al. 1993;Gorsevski et al. 2006).
Probabilistic methods such as LR (Liao et al. 1988), Bayesian analysis (Juang et al. 2001;Juang et al. 2002), Bayesian regression (Cetin et al. 2002;Cetin et al. 2004;Moss et al. 2006) and Artificial Neural Networks (Pal 2006;Samui et al. 2008;Xue and Yang 2014) have been recently applied in liquefaction studies at site scale to evaluate the factor of safety against liquefaction of a vertical soil column. These approaches, as well as commonly used deterministic methods, require a large amount of site-investigations in order to be applied over large areas. Indeed, the assessment of the factor of safety against liquefaction in a uniform soil layer in geoengineering practice is based on the ratio between cyclic resistance ratio (CRR) of the soil and the cyclic stress ratio (CSR) induced by design seismic acceleration (Youd and Idriss 2001;Seed et al. 2003;Cetin et al. 2004;Idriss and Boulanger 2008). The CRR is defined using in-site investigations such as cone, standard and Becker penetration tests or shear wave velocity determinations (Youd and Idriss 2001;Jafarian et al. 2010;Monaco 2014). Similarly, to compute the liquefaction potential index (LPI) of stratified deposits (Iwasaki et al. 1978), the factor of safety of each single layer must be assessed (Papathanassiou et al. 2015). Thus, over large areas, the application of deterministic or probabilistic methods to assess LPI is a demanding task in terms of the large number of in situ tests needed to assess the spatial variability of geotechnical and hydrogeological properties of the subsoil. Only multiple borings or soundings allow eventually to cope with under-sampled locations by using different interpolation methods (Baise et al. 2006;Holzer et al. 2006;Chen et al. 2008;Facciorusso et al. 2010;Pokhrel et al. 2012).
Considering that sectors of the Po Plain have not been investigated by seismic microzonation studies in all areas, the assessment and mapping of CSE susceptibility by means of spatial geostatistical approach is believed to be a well-justified option for contributing to the regional seismic risk planning and management.

Study area and coseismic surface effects of 2012
The study area extends for 1480 km 2 in the Po Plain, inside the Emilia Romagna Region (Northern Italy). It includes 19 municipalities within the provinces of Ferrara, Bologna and Modena (Figure 1). This area is economically well developed, with a highly specialized industrial sector, intensive agriculture and high population density. The hydrography of the area is defined by the rivers Po, Panaro, Reno and Secchia, flowing through the alluvial flood plain at altitude ranging from +32.3 m a.s.l to ¡0.2 m, with an average topographic slope lower than 3 .
In May 2012, the area was affected by a seismic sequence that peaked to Mw = 5.9. The main shocks of May 20 and May 29 induced hundreds of CSEs that were directly or indirectly related to liquefaction phenomena ( Figure 1).
The seismicity of this region is related to active tectonics in the buried front of the Apennines, which is characterized by faulted-folds (known as Ferrara Folds), forming an N-verging arc located right below the Plio-Quaternary marine and alluvial sediments of the Po Plain (Pieri and Groppi 1981). The structures responsible for the 2012 earthquake were the most external active Ferrara and Mirandola thrusts and folds (Bignami et al. 2012). The seismic sequence of 2012 occurred in the upper 10 km of the crust (Alessio et al. 2012).
The CSEs occurred during the 2012 seismic sequence have been surveyed by different teams of researchers and have been included into different open access maps and databases (Servizio Geologico Sismico e dei Suoli 2012; Emergeo Working Group 2013; ISPRA 2013). These reports include phenomena such as fissures, sand boils, ejections of sand and water from wells, bulging, lateral spread (Figure 2), as well as hydrological anomalies such as water table fluctuations and changing of water temperature or spill of wet sand from wells. Since Roman times, the natural and artificial diversions of river courses influence the geomorphology of the Po plain, shaping landforms particularly susceptible to CSE. The subsoil lithology, down to a depth of about 20-30 m, is characterized by poorly consolidated fluvial deposits dated from the Upper Pleistocene to the Holocene. More in detail, the uppermost tens of meters are essentially constituted by sandy and silty-sand deposits, which are spatially related with river channels, levees and crevasse splays, representing the paleo-drainage patterns over the floodplain Civico et al. 2015). Data from ARPA Emilia Romagna (i.e. the environmental agency of Emilia Romagna that monitors hundreds of wells in the area) indicate that the level of groundwater in the study area ranges from ¡3.5 m to maximum ¡16 m below the topographic surface.
Although recent technological developments have facilitated the application of spatial statistical modelling (van Westen et al. 1997;Carrara et al. 1999;Guzzetti et al. 1999;Guzzetti 2005), the implementation of such methods requires an important preliminary phase of data acquisition and information processing aimed to the proper selection of supporting evidences and of causal factors to be considered. In this research, susceptibility modelling with WofE and LR is based on the use of a subset of all the CSE occurred during the 2012 earthquake (as supporting evidence) and of a set of maps representing independent predisposing factors. The considered predisposing factors are static local terrain conditions, such as lithological, hydrogeological and geomorphological conditions that, according to literature, play the bigger role in determining CSE (Martelli 2012;Monaco 2014). Moreover, the absolute elevation of topographic surface and the distance to faults were also considered as predisposing factors. On the other hand, the 2012 earthquake seismic acceleration was not considered, since it is indeed a triggering parameter that, as such, should by no means be included in a susceptibility assessment. Another subset of CSE has then been used for model validation using Prediction Rate Curves (PRCs) and the Ratio of Effectiveness (RofE) (Chung and Fabbri 2003), as detailed in Section 3.4.
All thematic maps regarding supporting evidences and predisposing factors were generated using ArcGIS Desktop 10.3 software. The WofE and LR models were run in the freeware Arc-SDM Spatial Data Modeller extension (Sawatzky et al. 2004). Both CSE evidences (inventory dataset) and the predisposing factors maps (causal factors dataset) were converted from vector to raster structure with a cell size of 5 £ 5 m and projected in WGS 1984 UTM Zone 32 N.

Dataset: CSE inventory map
The liquefaction-related CSE inventory map ( Figure 1) was obtained by merging information from three main databases: Servizio Geologico Sismico e dei Suoli (2012) For modelling purposes, the conversion of 1910 point features into raster format resulted in only 1875 cells with a CSE evidences attribute. The CSE inventory map was subdivided in two sub-datasets: a modelling sub-dataset with 1306 CSE supporting evidences (70%) used to train WofE and LR models, and a validation sub-dataset with 569 CSE supporting evidences (30%) used to test and compare the predictive performances of the modelling results. The sub-datasets are the result of different tests of random partitions, which were continued until the percentage of the CSE sampled in the modelling sub-dataset provided a well-distributed representation of CSE occurrences.

Dataset: predisposing factors maps
The selected predisposing factors for CSE susceptibility modelling were the following: geomorphological units; lithological units; water table depth; topographic surface elevation; distance to faults and overthrusts (see also Table 1). The predisposing factors were classified in categorical layers and further reclassified into integer format to be included in the modelling process.
The geomorphological units were digitized from a geomorphological map of the Po Plain 'Carta Geomorfologica di Pianura' at 1:250,000 scale (Castiglioni et al. 1997). The geomorphological map used in this study is an adaptation of the previous map with more detailed data. From the original analogical map, the main fluvial geomorphological forms were digitalized in vector format and then converted to raster structure. The five geomorphology classes considered were: (1) crevasse splay; (2) levee well defined; (3) levee less defined; (4) abandoned riverbed; and (5) no fluvial landform (Figure 3(a)).
The lithological units were defined based on the official geologic map 'Carta Geologica di pianura dell'Emilia-Romagna 1:250,000' (Servizio Geologico, Sismico e dei Suoli 1999). Specifically, eight lithological units were considered (Figure 3(b); Table 1): (1) medium fine sand, in beds tens of cm thick; (2) sandy silt; (3) silty clay; (4) silt and clayey silt; (5) clayey silt and sandy silt; (6) medium to fine sand; (7) medium fine sand, silt and silty clay; (8) medium and coarse sand. Lithological units 6 and 7 belong to deltaic and coastal deposits, whereas unit 9 corresponds to the Po riverbed that is part of the study area. All remaining lithological units are associated with the alluvial plain deposits. Groundwater depth in the phreatic and upper confined aquifers was mapped using a dataset provided by the Regional Agency for Prevention, Environment and Energy of Emilia-Romagna region (ARPAE) based on a network of hundreds monitored wells. Data regarding the year 2012 were rasterized by interpolation of 1 m isolines. Considering that the maximum groundwater depth associated to the possible occurrence of CSE is about 16 m in these phreatic and upper confined aquifers Monaco 2014), the water table depth was divided into seven classes: (1) 3.5-4.5 m; (2) 4.5-5.5 m; (3) 5.5-6.5 m; (4) 6.5-7.5 m; (5) 7.5-10.5 m; (6) 10-13.5 m; (7) 13.5-16.5 m (Figure 3(c)).

Formulations of the models
The WofE is a simple statistical method based on Bayes' theorem (Lee 1989;Agterberg et al 1993;Bonham-Carter 1994;Denison et al. 2002). The method is applicable when sufficient data are available to estimate the relative importance of evidential themes via statistical means (Bonham-Carter 1994) computing prior (unconditional) probability and posterior (conditional) probability. The method computes the conditional probability that an event (in this case CSE) does or does not belong to a given set of overplayed causal factor maps. If a casual factor is B, the classes in it are Bi and supporting CSE evidence is (s), and using the Bayes' theorem, it is possible to calculate the conditional probability for each class Bi as follows:  Table 1. where P(Bijs) is the conditional probability to have Bi given s, P(s) is the prior probability to find s within the study area (SA), and PBi is the prior probability to find the class Bi within the study area. Then, the conditional probability to have s when the class Bi is not present is where P(Bi^js) is the conditional probability not to have the class Bi given s, P(s) is the prior probability to find s within the study area AS, and P(Bi^) is the prior probability not to find the class Bi within the study area.
In the WofE, the weight of every class of causal factors is calculated by a combination of positive and negative correlation values defined with Bayes method (W+ and W¡) and by the difference between them (contrast 'C'). This method assigns positive and negative weights (Wi+ and Wi¡) to each pixel of the predisposing factors maps.
The difference between the two weights is known as the weight contrast C. The magnitude of the contrast reflects the overall spatial association between the predictable variable and the CSE.
A value of contrast equal to zero indicates that the considered class of causal factor is not significant for the analysis. A positive (negative) contrast indicates a positive (negative) spatial correlation between factor class and CSE occurrence.
The susceptibility map is obtained by a linear combination of weights of different factor classes. Then, the value of posterior probability Pp(s) is given by the following formula: The LR is a multivariate statistical method suitable for estimating the probability of a dichotomous dependent variable, such as the occurrence or absence of the phenomenon, from its relations with independent variables (Cox 1958;Agterberg et al. 1993). The method considers several physical parameters (variables) that may affect the probability of an event and yields coefficients for each variable based on sample data. The algorithm of logistic regression applies maximum likelihood estimation after transforming the dependent variable into a logit variable (natural log of the odds of the variable occurring or not). In this way, LR estimates the probability of a certain event occurring. The resulting probability values express the degree of CSE susceptibility, and they range between 0 and 1. LR can be applied when the variables show conditional dependence, contrary to the WofE method. It can also be performed using categorical or continuous variables as well if they are not normally distributed (Hosmer and Lemeshow 1989). The correlation between the CSE evidences in a terrain unit area and the casual factors conjunction is expressed by the following equation: where P is the probability that a unit area falls into CSE or no CSE liquefaction condition. LR fits a special s-shaped curve by taking a linear regression, that may produce any y-value between ¡1 and +1, and transforming it with the function that produces a probability (p probability) between 0 (as y approaches minus infinity) and 1 (as y approaches plus infinity). Z is the linear combination of: where b0, b1, … are the regression parameters; x1(r), x2(r), …, xm(r) are the independent variables for each cartographic unit; e is the error associated to the curvilinear approximation of the model. In the logistic function (Equation (9)), the input is z and the output is P. The variable z represents the exposure to some set of independent variables (xm), whereas P represents the estimated probability of CSE, given that set of explanatory variables. It must be noted that the LR approach does not require, or assume, linear dependencies between the dependent term and the variables considered (independent set of variables representing causal factors). An exponential function is involved. Coefficients b are estimated through the maximum likelihood criterion and correspond to the estimation of the more likely unknown factors. In LR, the coefficients b express the relative importance of the classes. A positive b coefficient indicates that the class contributes to susceptibility, whereas a negative b coefficient class reduces CSE susceptibility. Coefficients b close to zero indicate little relation between the class and the response variable (e.g. CSE) and coefficients equal to zero indicate the class is not influencing the susceptibility.

Models application and validation procedures
The conditional independence existing among predisposing factors is a prerequisite of the WofE and LR methods to avoid overestimations of the posterior probability (Blahut et al. 2010;Pereira et al. 2012). In this work, we assume that auto-correlation among variables can exist. For that reason, final results are not considered in probabilistic terms but only to rank susceptibility within the study area. The susceptibility classes are defined following the same principle.
The predictive capability of the outputs of the WofE and LR models was assessed by using the PRC and Success Rate Curve (SRC) (Chung and Fabbri 2003;Schmitt 2010) and by computing their area under the curve (AUC) (Begueria 2006;Zêzere et al. 2007). Another procedure herein used to estimate the real predictive capability of each class of susceptibility and then for model comparison was the RofE (Chung and Fabbri 2003) of CSE susceptibility classes. Chung and Fabbri (2003) pointed out that a RofE equal to 1 indicates a random and unsupportive class prediction. Good class power predictions correspond to RofE higher than 3 (less than 0.2) for high (low) susceptibility classes. Hence, the numeric representations can be assumed as quantitative measures of the predictive power of the class when the probability function is used to assess the susceptibility of future CSE.
In order to compare the CSE susceptibility maps obtained by the two models, the susceptibility maps have been classified using the same number of classes, defined by ranking results in decreasing order of susceptibility and by fixing the limits between classes by using the same values, for all models of cumulated spatial extent. The limited scores obtained for the posterior probability provided by LR have forced us to subdivide the study area into only four classes, where the first 30% was discriminated in three classes representing each 10% of the total study area. The fourth class occupies the 70% of the target area and corresponds to the low and very low susceptibility.
The degree of spatial correlation and the relative degree of class dispersion between CSE susceptibility maps obtained with WofE and LR was estimated on the basis of the Cohen's kappa coefficient (Cohen 1960;Hoehler 2000).
Finally, to highlight the extent to which CSE susceptibility should be considered more certain, the output of the two models were merged so to obtain a final synthetic map of the more susceptible zones that highlights which areas are considered at a high level of susceptibility by one or both models.

Distribution of CSE with respect to predisposing factors
The distribution of the CSE in the different classes of each predisposing factor is represented in Table 1. Regarding geomorphological features, it should be highlighted that fluvial landforms covering 27% of the study area include 60% of all the CSEs (Figure 3(a) and Table 1). As can been seen, the majority of CSEs were distributed along the main ancient fluvial landforms. This finding is consistent with the fact that levees, crevasse splays or abandoned riverbed or land reclamation areas generally act as preferential location for the CSE occurrence .
The sub-surface lithology seems also to play a discriminant role, as almost 60% of CSE are associated to the medium fine sand of alluvial deposits (1) and 32.7% to the silty clay lithology class (3) (Figure 3(b) and Table 1). The CSE distribution seams also directly related to groundwater depth. In the study area, 35% of CSE are associated to groundwater depth of 3.5-4.5 meters and another 53% of CSE is located in areas with groundwater depth between 5.5 and 7.5 m (Figure 3(c) and Table 1). This is consistent with the observation of Monaco (2014), which indicates 15 m as the maximum groundwater depth associated to liquefaction phenomena. Despite the almost flat morphology of the Po Plain, in order to verify any possible influence of elevation on CSE occurrence, the DEM was also included in this evaluation. The elevation ranges of 12-15 m a.s.l. has the highest percentage of CSE (36.8%) (Figure 3(d) and Table 1). This elevation is generally associated to the levee of the ancient riverbeds and fluvial depositions forms. For what concerns the distance to faults, the 55.5% of CSE occurs in the class 0-500 m, i.e. very close to tectonic structures. The remaining 44.5% is distributed in progressively lower percentages as distance to tectonic structures increases (Figure 3(e) and Table 1).

CSE susceptibility using WofE and LR models
The application of the WofE method resulted in the positive and negative weights (W+ and W¡) and, consequently, in the weights contrasts (C) listed in Table 2. Regarding geomorphological features, the higher positive values of C are obtained with respect to well-or less well-defined levees and crevasses splays. Concerning lithology, the higher positive values of C are obtained for medium fine sand, whereas significant negative C is obtained for sandy silt and the Po riverbed deposits. As for groundwater depth, the higher positive values of C correspond to the range of depth 6.5-7.5 m and the maximum negative C to the class 10-13.5 m below the ground. Other classes are poorly characterized by the contrast. Elevations class from 12 to 15 m shows the higher positive C value. The distance to faults shows the higher positive C value in the class from 0 to 500 m, and the C becomes significantly negative for classes representing distances higher than 2 km.
The application of the LR method resulted in the b coefficients summarized in Table 2. Since the b coefficients indicate the importance of each variable in the post-probability calculation, it can be stated that the lithology, again, is the most important conditioning factor. In particular, the lithology class 'medium fine sand' shows the highest b coefficient. Positive values are obtained also for other lithological units: sandy silt, silty clay, medium and coarse sand and Po riverbed deposits. Depth to groundwater presents positive regression coefficients for all the classes. On the other hand, the lowest values of b coefficients are associated to the distance to faults and to the DEM. Nevertheless, the significant negative values for classes of distances to fault higher than 2 km indicate that these classes are inversely correlated to the CSE occurrences.
The application of the WofE method to the entire dataset of predisposing factors returned posterior probability values ranging from 0 to 0.507, whereas the application of LR generated posterior probability values ranging from 0 to a maximum of 0.252.
In order to retrieve classified susceptibility maps, the WofE and LR probability values have been divided into four classes by ordering probability values in decreasing order. The higher susceptibility classes ('very high,' 'high' and 'moderate') are defined at fixed 10% intervals of the study area. The 'low to very low' class is assigned to the low probability values obtained in the remaining 70% of the  Figure 4), the medium, high and very high susceptibility zones are concentrated along levees or areas characterized by medium to fine sand, in the districts of Bondeno to Mirandola and Concordia, as well as San Felice and Sant'Agostino ( Figure 4).

CSE susceptibility models validation and comparison
The two statistical models were validated and compared through the PRC computed using the validation subset (that includes 569 CSE). Both methods show very similar results, with AUC of the PRC of 0.83 for WofE and 0.85 for LR, and slightly lower values for the AUC of the SRC of both methods ( Figure 5). According to Guzzetti et al. (2006), when the overall model returns an AUC greater than 0.8 (i.e. 80%), the model can be considered an acceptable predictor. Looking at the trend and benchmark values of the PRC and SRC ( Figure 5), both models predict about 60% of validation CSE within 10% of area, 70% of validation CSE in 20% of the study area and 90% of validation CSE in 50% of the study area. In both models, the values of RofE corresponding to the very high susceptibility class is higher than 6, thus much higher than for the other susceptibility classes that, indeed, remain lower than 1 (Table 3). This, also, is an indicator of good quality of the prediction. The degree of spatial correspondence between CSE susceptibility maps was assessed using the Cohen's kappa coefficient (K) and on the total spatial agreement. Despite the moderate K value (K = 0.68), the total spatial agreement within the study area reaches 85.67% (Table 4). Specifically, the very-high, high, moderate and low to very low susceptibility classes in the WofE and LR models returned a spatial agreement in the order of 80%, 53%, 50% and 96%, respectively. This highlights some differences in the spatial distribution of susceptibility classes between the two models, especially in the high and moderate classes. Consequently, an attempt has been made to create a synthesis map that shows which areas are classified as being at very high, high and moderate susceptibility by both models. This resulted in a final map that highlights the zones indicated as susceptible by both models should be considered more probably affected by future CSE occurrences ( Figure 6). In practice, the area classified from very high to moderate susceptibility by both the models covers 24.6% of the study area and includes 73% of the inventoried CSE of May 2012.

Discussion
The aim of the study was to compare bivariate and multivariate analyses with two different datadriven methods (WofE and LR) in order to evaluate liquefaction-based CSE susceptibility in the central part of the Po Plain. Here, on May 2012, the two main shocks, of a long seismic sequence, induced widespread CSE near the epicentral area. The outcomes obtained by the parallel application of WofE and LR data-driven methods raise some points of discussion. Both methods have provided acceptable results, expressed by PRCs, obtained using an independent validation dataset, showing an AUC above 80%.
Both models indicated that the spatial distribution of susceptible zones is mostly controlled by fluvial landforms as crevasse splays, fluvial levees, and ancient riverbeds. This is in agreement with previous studies that considered the presence of paleo riverbeds, out-of channels fans among the main controlling factor for the occurrence of CSE during the 2012 seismic sequence (Papathanassiou et al. 2015). The medium to fine sandy lithology associated to these facies beside the shallow water table depth represents the other main casual factors which are outlined in the distribution of susceptibility areas. The elevation of ground surface provides a positive contribution to susceptibility prediction just in one class of elevation that for the southern part of the study area appears associated mainly to the levees of fluvial geomorphological landforms. As shown in the output maps (Figure 4), the susceptible zones cover the central part of the study area, displaying two wide extreme areas (east and west zones) of low and very low susceptibility for all the maps produced. On the other hand, the central classes of susceptibility seemed to be more spatially scattered, showing also a higher degree of discrepancy between the models' classification.

Conclusion
The existing body of research and the literature produced after the May 2012 earthquakes focused to assess the liquefaction potential at municipality level with deterministic methods. CSE-liquefaction susceptibility assessments evaluated with geostatistical approach remains an open research field rarely tested.
During the last few decades, both deterministic and probabilistic conceptual understanding and procedures have remained strictly centred on liquefaction potential evaluations performed in situ. Starting from site evaluations, the susceptibility assessment over large areas is generally based on interpolations methods, which contain intrinsic not liquefaction phenomena-related generalization processes. In this sense, this work demonstrates that CSE susceptibility can be assessed at the regional scale based on the statistical relationship between CSE evidences and a dataset of terrain variables assumed as CSE predisposing factors with good accuracy using bivariate and multivariate geostatistical data-driven methods.
From a methodological perspective, this study has pretended to propose a contribution in the branch of CSE susceptibility evaluations at regional level. Considering the results obtained and the quality of the susceptibility maps produced within this study, we conclude that WofE and LR can be applied as assertive methods for susceptibility evaluations over large areas as the first level of liquefaction analysis. However, the predictive capacity of CSE susceptibility models is strongly dependent on the completeness and accuracy of the CSE inventory. Recently, the use of differential synthetic aperture radar (DInSAR) coherence data has been proved to be effective to map the surface liquefactions effects in Japan (Chini et al. 2013) and in Italy (Chini et al. 2015), which can be used to improve the susceptibility models in future work. Besides the validated accuracy of the obtained results, geostatistical data-driven methods demonstrated to be reliable, particularly for their low operational cost and promptness of procedures. One of the advantages of the CSE susceptibility maps produced with the statistical approach deals with the capability to assess the susceptibility to CSE for all the study area, including also those areas that are not taken into consideration by the official microzonation cartography (areas outside urbanization planning). In addition, throughout the susceptibility maps, the applied approach allows to promptly fulfil a synthetically overview of the phenomena at regional scale. Furthermore, the outputs consent to spatially differentiate the phenomena into different degrees of susceptibility within the municipality boundaries. Considering the degree maps detail, the susceptibility delineated resulted to be more precise than the first level of seismic microzonation studies (e.g. PTCP). Although by a conservative point of view these maps cannot be directly applied to land-use restrictions, they represent a useful cartographic reference during emergency or in planning long-term risk reduction strategies and should be considered as a part of the background knowledge documentation within the official seismic risk cartography.
These supporting reasons suggest testing the approach presented in this study in other seismic areas beyond Emilia, characterized by similar geological and fluvial-dominated environments and potential liquefaction susceptibility proneness.