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