Settlements and angular distortions of shallow foundations on liquefiable soil

The extensive injury caused to buildings by liquefaction during past earthquakes, with the uneven spatial distribution of damage, raises the need for rapid predictive tools, applicable at a large scale, that comprehensively account for the properties of earthquake, subsoil and structure. A method is herein proposed to quantify the angular distortion of framed low-rise buildings based on a simple characterization of the above factors. The analysis moves from past literature criteria introduced to quantify the vulnerability of buildings under static conditions and extends their applicability to liquefaction assessment integrating parametric two-dimensional numerical analyses with recent literature predictive formulas and machine learning inference. The numerical calculation, performed for variable stratigraphic and mechanical characteristics of the subsoil, ground motion and equivalent flexural stiffness of the foundation, quantifies the role of each factor on the absolute settlement and angular distortion. Then the dependency on the different factors of the angular distortion is inferred with an artificial neural network, grouping parameters to limit the number of input variables and express results with charts that make prediction more accessible.


Introduction
Seismic liquefaction continuously stimulates the interest of the geotechnical community in promoting experimental, theoretical studies and field observations that clarify the factors ruling susceptibility of soil, triggering (e.g., Ishihara 1996;Youd and Idriss 2001;Seed et al. 2003;Bray and Sancio 2006;Boulanger and Idriss 2014) and chained mechanisms (Iwasaki et al. 1978;van Ballegooy et al. 2014;Chiaradonna et al. 2020;Cubrinovski and van Ballegooy 2017). For buildings and infrastructures, studies are oriented at predicting damage (Bird et al. 2005;Bray and Macedo 2017;Karamitros et al. 2013;Bullock et al. 2018;Castiglia et al. 2020; Baris et al. 2021) or at conceiving mitigative solutions (Flora et al. 2021;Salvatore et al. 2020), motivated by the concern of stakeholders for the huge economic and social impact of liquefaction recorded in past seismic events. An overall estimate is provided by Daniell and Vervaeck (2012) who disaggregated primary (shaking) and secondary causes (tsunami, fire, landslides, liquefaction, fault rupture, and other type losses) over a global record of nearly seven thousand earthquakes from 1900 to 2012. These authors found liquefaction responsible for about 2.2% of the direct economic losses, globally estimated at 2.24 trillion US dollars, this fraction becoming 3.6% when considering total losses, i.e. direct plus indirect damage. This relatively small percentage could misleadingly drive to underestimate the relevance of liquefaction, but the numerous examples that occurred in urbanized systems prove that the physical damage together with the prolonged impracticability of buildings and infrastructures can be overwhelming, as it undermines the recovery of normal living conditions and may lead to the complete abandonment of the place (Macaulay et al. 2009;CSAPEISLA, 2016). Literature reports plenty of examples of destructive effects caused by liquefaction on urban areas (San Francisco 1960-Youd andHoose 1978;Kobe 1995-Chung 1996Kocaeli 1999-Cetin et al. 2002Christchurch 2010-2011-Cubrinovski et al. 2011Urayasu 2011-Yasuda et al. 2012; Baris et al. 2021;Emilia-Romagna 2012-Fioravante et al. 2013. The soil-foundation-structure interaction strongly rules this phenomenon . Fully coupled analyses of liquefaction and shaking are rather complex and require a very accurate definition of the input variables, difficult to achieve when performing serial analyses like in large-scale assessment. Bird et al. (2006) propose a simpler solution that decouples and recombines the effects of shaking and liquefaction for a prescribed event. The first compulsory step is the definition of a demand variable, i.e. the most expressive quantity that characterizes physical damage. Foundation movement is classically subdivided into the following components: mean settlement, rigid rotation and angular distortion. Rotation is of paramount importance for high-rise buildings and towers, prEN 1997-1 (2004), both in static and seismic cases. Poulos (2001) converge on identifying the angular distortion as the paramount demand variable for low-rise buildings (1-3 stories, Fotopoulou et al. 2018), fixing thresholds equal to 1/500, 1/300 or 1/150 for respectively light, visible cracking, and structural damage of framed buildings in the static case. A similar assumption is made by previous studies that considered the effect caused by the self-weight of buildings (e.g. Skempton and MacDonald 1956;Bjerrum 1963;Meyerhof 1953;Burland and Wroth 1974;Burland et al. 1977;Wahls 1981). Boscardin and Cording (1989) studied the effects of excavation and defined the damage level as a function of angular distortion and horizontal tensile strain, the two quantities related to each other for the different sources of movement. This step is simplified by Grant et al. (1974) that propose to statistically estimate angular distortion as the upper bound from a relation with absolute settlements measured on buildings founded on cohesive or cohesionless soils with rafts, strips or isolated footings.
For liquefaction assessment of low-rise buildings, Bird et al. (2006) distinguish the case of rigid from flexible foundation, with absolute settlements being the principal damage factor for the former, differential settlement for the latter case. In this circumstance, differential settlements induced by liquefaction on framed buildings cause a drift of columns that the authors propose to cumulate with that produced by shaking. Fotopoulou et al. (2018) also adopted the differential settlement as a demand variable in their probabilistic definition of vulnerability for low-grade structures. In general, the structural damage induced by foundation movements on a building depends on the stiffness and fragility of the structurefoundation system, these factors are connected with the typology, extension, and height of the building. A continuous transition, rather than a net separation between rigid and flexible structures, would better describe the variety of possible structural typologies.
For liquefaction, several methods are proposed to predict absolute settlements. The state-of-the-art practice for this evaluation largely relies on empirical procedures developed to estimate post-liquefaction, one-dimensional consolidation settlement in free-field conditions (e.g., Tokimatsu and Seed 1987;Ishihara and Yoshimine 1992). According to Ishii and Tokimatsu (1988), this assumption can be reasonably accepted only if the width of the foundation is at least twice or three times larger than the thickness of the liquefiable soil layer. However, the main limitation of these empirical procedures is that none of them considers the soil-structure interaction and the resulting complex mechanisms, for example, the SSI-induced building ratcheting during earthquake loading (Dashti and Bray 2013). Based on the results of numerical analyses and attributing liquefaction-induced settlements to the seismic excitation characteristics and the post-shaking degraded static factor of safety, Karamitros et al. (2013) provide a simplified analytical formula for the estimation of absolute settlement of strip and rectangular footings with a clay crust. Such settlement is associated with a ''sliding-block'' type of punching failure through the clay crust and within the liquefied sand layer. Bray and Macedo (2017) performed a large number of parametric numerical analyses and proposed to express the total settlement as a sum of three contributions respectively induced by shear, volume deformation and sand ejecta. In particular, the shear-induced rate is related to several properties, including the unitary contact pressure on the foundation, the thickness of the liquefiable layer and the lower planimetric dimension of the building footprint and the cumulative absolute velocity (Campbell and Bozorgnia 2011). Performing rich and various parametric numerical, fully coupled three-dimensional analyses of the soil-structure interaction, Bullock et al. (2018) define a relation to predicting the statistical distribution of settlements for shallow-founded structures on liquefiable soil induced by volumetric and distortional strains. This formula has the advantage of capturing the role of most soil, ground motion and building properties.
The present paper summarizes the above concepts proposing a rapid calculation of the angular distortion induced by liquefaction on shallow-founded low-rise buildings, applicable for the preliminary screening of risk at a large scale, e.g. cities, and districts. To this aim, a large number of two-dimensional coupled numerical analyses are performed, parametrically varying the stratigraphic and mechanical conditions of the subsoil, the ground motion characteristics and the equivalent structural properties of the building. The analysis leads to inferring a relation with the relevant properties using artificial neural networks.

Definition of the model
The implemented numerical model simulates with a Finite Difference code (FLAC v8, Itasca, 2016) the two-dimensional layout depicted in Fig. 1 developed over a width of 40 m and a depth of 20 m. Luque andBray (2015, 2017) showed that the primary aspects of the dynamic response of a 3D system in terms of liquefaction-induced building settlements can be captured by 2D analyses when tributary masses and stiffnesses are given to the buildings. The three-layer subsoil model, that includes an upper low-permeability cap (Layer #1), an intermediate liquefiable layer (Layer #2) and a stiffer base (Layer #3), has been chosen as one of the most frequent subsoil conditions. The predominant role is played by the intermediate liquefiable layer. Tests are available in the literature (e.g. Millen et al. 2020) to assess whether real conditions conform to this scheme characterized by a single liquefiable layer or alternative models including multiple liquefiable layers are needed. The calculation mesh consists of 13,980 rectangular elements of 0.8 m in width and variable heights (0.5 m for the above and below layers, 0.4 for the liquefiable layer). These dimensions were chosen after the suggestion of Kuhleimeyer and Lysmer (1973), who found that the propagation of seismic waves in continuum media can be simulated with sufficient accuracy if the element's dimension is smaller than 1/10 of the minimum propagating wavelength.
The stress-strain response of the liquefiable soil has been simulated with PM4Sand Version 3.1 (Boulanger and Ziotopoulou 2017), thanks to its capability to capture the cyclic behaviour of saturated sand. The stress-strain response of the shallow and deepest layers is simulated using Mohr-Coulomb hysteretic model, considering the stress state induced into the subsoil by the above building. In principle, upper and lower layers affect the seismic excitation of the whole system and concur to determine the settlements of foundation. Thus, more refined constitutive models would be preferable for these layers too. For instance, the pore pressure build-up in the upper crust could reduce shear strength and eventually increase settlements. To this aim, an interesting analysis is provided by Dahl (2011), who compared the cyclic undrained shear strength t cyc of materials with different plasticity to the static undrained shear strength s u . Their study shows that t cyc is lower than s u with differences in the order of 30 or 40% (with respectively 10 or 30 cycles) for low plasticity soils (PI < 10%). For higher plasticity materials, differences are less relevant, being in the order of 10 to 20% for respectively 10 or 30 cycles. Therefore, the use of more complex models (e.g. PM4Silt-Boulanger and Ziotopoulou 2018) would be justified for low plasticity soil, while would introduce an unnecessary complexity for more plastic soils. In the present study, crust broadly indicates a category of plastic soils. Furthermore, the adopted simplification, made in the spirit of a preliminary calculation, has the advantage that the constitutive parameters of Mohr-Coulomb model can be determined with routine in-situ or laboratory tests. Nevertheless, the variability of seismic excitation in the intermediate layer is investigated considering largely different seismic inputs in the parametric analysis. The constitutive parameters for each model (see Table 1) have been set from the back analysis of a case study in Terre del Reno (Italy) (Fioravante et al. 2013;Sinatra and Foti 2015;Facciorusso et al. 2016). In this scheme, the upper crust is made of silty-clayey soil (PI = 15 ÷ 20%), while the lower base is made of relatively stiff clays. The availability of data on a such well-documented case study enabled the validation of the adopted mechanical schematization versus observation (see "Appendix 1"). The calibration on this specific case study does not affect generality, as the role of the most relevant soil properties has been explored by systematically modifying them in a parametric analysis.
Starting from the natural unit weights and the void ratios reported for the different strata by Fioravante et al. (2013), the solid phase and dry unit weights are assigned in Table 1. In the same table, permeability for each stratum is derived from the work of Sinatra and Foti (2015) based on the analysis of CPTU tests. The small strain elastic stiffness parameters are inferred from the shear wave propagation velocity by computing the small strain shear modulus before and computing the bulk modulus by imposing a Poisson's ratio equal to 0.3.
Stiffness degradation and damping for the lower base and crust have been simulated through a cubic equation with two parameters c 0 and c 1 calibrated by fitting the experimental G/G 0 and D vs. γ curve provided by Fioravante et al. (2013) and Sinatra and Foti (2015) (Fig. 2). In this particular calculation, cyclic simple shear tests have been simulated with the numerical code and c 0 and c 1 have been fixed with a trial-and-error procedure looking for the best match between numerical and experimental results.
The strength parameters for the Mohr-Coulomb model, namely friction angle (ϕ'), cohesion (c') and undrained shear strength (c u ) are extracted from the work of Sinatra and Foti (2015) that provides them from the analysis of laboratory and CPT tests. The effective strength values and the large strain stiffness parameters for the cohesive soils are taken from the literature (Itasca Consulting Group, Inc. 2016).
The PM4Sand Version 3.1 (Boulanger and Ziotopoulou 2017) adopted for the liquefiable layer represents an evolution of the Dafalias and Manzari (2004) formulation, being a stress-ratio controlled, critical state compatible, bounding surface plasticity model developed primarily for earthquake engineering applications (Boulanger andZiotopoulou 2013, 2015). In PM4Sand the Jefferies state variable = e − e cs is modified as proposed by Boulanger (2003): where D R and D R,CS are the relative densities at the current and critical states, respectively. Deformation is regulated by the stress ratio and develops according to multiple bounding surfaces. In particular, the model incorporates the bounding, dilation and critical surfaces of Dafalias and Manzari (2004) but removes the Lode angle dependency on friction angle. Thus the bounding (M b ) and dilation (M d ) ratios are related to the critical stress (M) ratio by the following simpler expressions: where n b and n d are calibration parameters and M is computed as a function of the critical state friction angle.
The PM4Sand model has been calibrated to capture the cyclic undrained behaviour of the liquefiable sandy layer seen by Facciorusso et al. (2016). In particular, four triaxial undrained cyclic tests have been numerically simulated with the adopted code (FLAC v8, Itasca, 2016) matching the liquefaction resistance curve (CSR vs. N liq ). The calibrated parameters are listed in Table 1 while the comparison of the numerical and experimental results is shown in Fig. 3.
The superstructure is modelled with an equivalent beam characterized by a flexural stiffness (EI) and a contact pressure (q). The EI modulus summarizes the flexural stiffness of Fig. 2 Shear modulus degradation and damping curves numerically and experimentally (Fioravante et al. 2013;Sinatra and Foti 2015) obtained for the different subsoil layers the foundation-superstructure complex but, for the sake of simplicity, this modulus can be conservatively computed for the sole foundation. The contact pressure q summarizes the contribution of all building floors and can be computed by multiplying the number of floors times the overall unit load, i.e. the load per unit area including the self-weight of structural and non-structural elements and the accidental loads.
In "Appendix 1", the calculation has been performed for the considered case study, on one side showing a typical output of the numerical analysis, on the other side validating the calculation scheme with reference to a case where seismic input, subsoil and building conditions are known together with the surveyed settlements induced by liquefaction.

Parametric study
The physical, mechanical and geometrical factors varied in the parametric analysis have been chosen following Karimi et al. (2018). They observed that the mean permanent settlements of buildings are affected by contact pressure, seismic input, thickness, relative density and depth of the liquefiable layer, but also by the presence of a low-permeability cap. In this calculation, the parameters for the stiffer base (Layer #3) have been set equal to the reference case (Table 1), while the crust's strength has been varied parametrically giving values of undrained shear strength equal to 25, 50 and 100 kPa. This choice has been made featuring the three typical situations of respectively normally to slightly, moderately, and highly consolidated soils. Accordingly, E u /s u ratios equal to 550, 300 and 150 have been respectively assigned following the suggestion of Duncan and Buchignani (1976). In any case, the effects on seismic excitation have been considered by assigning largely different inputs, i.e. varying Arias Intensity and PGV.
The complete list of varied parameters is summarized in Table 2, positioning the water table at the ground level for all calculations. The structure-foundation system has been modelled with an equivalent continuous foundation system characterized by width (B), flexural stiffness modulus (EI) and embedment depth (H d ). The choice of continuous foundation is considered to simulate strip footings and mat foundations, which are frequently used for low-rise buildings on relatively weak soils. The assumed EI values (Table 1) correspond to reinforced concrete beams of different flexural stiffness. The maximum value (EI = 260MN/m) corresponds to a beam of rectangular Sect. 1.0 × 0.5 m (width x height), spanned with 4 m distance. Finally, the embedment depth has been varied between 1 and 3 m, which are typical values for shallow foundations (see Table 2). In addition, a set of six numerical analyses has been performed to evaluate the relevance of the structure's inertial mass, simulating height/width ratio (H/B) ranging between 0.5 and 1.5. Finally, considering the fundamental role played by the earthquake magnitude, four waveforms have been applied in the analysis, extracting them from the PEER Strong Ground Motion Databases. The velocity time history of these events, chosen thanks to their largely different Arias intensity (Table 3), have been scaled by three amplitude factors, respectively 0.7, 1.0 and 1.6 in order to explore the influence of a wider range of ground motions. This choice allows examining the influence of the intensity of the ground motion on liquefaction-induced foundation movements. Thus, a total number of twelve time-histories were assigned. Combining the seismic input with the parameters reported in Table 1, a total number of 353 analyses have been carried out. The sample calculation performed in "Appendix 1" shows that settlements start to develop for maximum r u values larger than 0.6 ÷ 0.8 in the liquefiable material (Layer #2). This observation, coupled with the summary of results reported in "Appendix 2" that shows maximum r u >0.8 for all cases, leads to conclude that the computed settlements are in all cases induced by liquefaction.
The typical output of calculation consists of the displacements profile below the foundation (Fig. 4) from which the following characteristic variables are extracted: • maximum, minimum and mean settlement: w Max, w min and w av ; • angular distortion: β; • horizontal displacement: S hi; • horizontal deformation: e h .

Sensitivity study
A sensitivity study has been initially performed to understand the relative influence of the varied parameters on the kinematic variables defined in Fig. 4. Firstly, the relation between mean and maximum settlements (w av and w MAX ) has been investigated. Their equivalence is needed for the following analyses where their outputs are alternatively related to the angular distortion of the building-foundation system. Figure 5 shows a proportionality between these two variables, being the ratio w av /w max equal on average to 0.84, with minimum and maximum values equal to respectively 0.73 and 0.98. The results of the calculation are then summarized by looking at the dependency of maximum, differential settlement and angular distortion (respectively w MAX , δ = w MAXw MIN and β defined in Fig. 4) on the most relevant parameters varied in the analysis. Noting that, within the investigated range (1 ÷ 3 m), embedment depth H d plays a lesser influence on the absolute and differential settlements, the role of the upper crust (Figs. 6, 7), liquefiable layer (Figs. 8,9) and structural characteristics (Figs. 10, 11, 12) is separately examined. parameters on all the considered components of the foundation movement, attenuation rates being more remarkable for the stronger seismic events. Settlements reduce almost linearly within the considered thickness range (H c ≤6 m) for the higher seismic intensities (f = 1.0 and 1.6), while reduction is smoother for the lower intensity earthquake (f = 0.7) (Fig. 6). Attenuation rate is very high for undrained shear strength s u increasing up to 50 kPa, then drops progressively for increasing s u (up to 100 kPa and more) (Fig. 7). In summary, despite preventing the excess pore pressure exhaust, the impervious crust forms a bridge below the foundation that contributes with its strength to limit the buildings settlements. Figures 8 and 9 point out the influence of the liquefiable layer. Here settlements and deformation increase rather continuously with thickness (H L ), being rates dependent on the earthquake intensity (Fig. 8). On the contrary, the variation determined by soil density is sharper. Effects are particularly critical for the lowest considered density (D r = 20%) with w MAX reaching values as high as 2.5 m (for f = 1.6), then sharply reduce for D r = 40% and D r = 60%. The latter result can be explained by the higher stiffness and lower tendency to contract seen on denser cohesionless materials (e.g. Modoni et al. 2000), while the coupled dependency on thickness and relative density justifies the assumptions made in the definition of liquefaction severity indexes (e.g. Iwasaki et al. 1978;van Ballegooy et al. 2014;Paolella et al. 2022).
Finally, Figs. 10, 11 and 12 summarize the effect of the building-foundation system. As expected, higher loads produce larger absolute and differential settlements (Fig. 10). Foundations width B produces a reduction on the absolute settlements throughout the considered range (Fig. 11a), less evident on the differential component (Fig. 11b, c). This effect, also seen by Bray and Macedo (2016), is somehow surprising if compared with the dependency of settlements induced by static loading. It can be explained considering the relevance of the shear stresses initially activated by the foundation loads in the subsoil on the triggering of liquefaction. This effect becomes progressively marginal with the foundation width, as this solution tends to the one-dimensional scheme characterized by lower shear stresses and higher mean effective stresses. Finally, flexural stiffness EI produces a continuous reduction in the differential settlements and angular distortion (Fig. 12b, c) but has negligible effects on the absolute settlements (Fig. 12a). As will be shown later, the above sensitivity analysis and the observed dependencies have been used to infer functional relations among foundation settlements, distortion and characteristic variables grouped in dimensionless parameters.

Role of building inertia
In the present study, the structure-foundation system has been simplified and modelled with an equivalent plate considering this case representative of low-rise buildings. This assumption has been made considering the observation of Karimi et al. (2018) on buildings of relatively limited height (H/B ≤ 3), who found a negligible role on the foundation settlements played by the structure's inertial mass and height/width ratio. For the sake of completeness, the consistency of this assumption has been herein verified by performing numerical analyses with B fixed equal to 10 m and varying the height/width ratio of the building (H/B) between 0.5 and 1.5. In this calculation, buildings have been simulated with equivalent linear elastic-perfectly plastic frames with pillars of 0.5 × 0.5 m section and beams of 0.3 × 0.5 m (width x height) section set at an inter-storey height equal to 2.5 m. A self-weight has been assigned to the beams equivalent to the unit load given in the simplified analysis. The effect of simplification is thus presented in Fig. 10, showing the ratio (ρ) between the foundation movements computed with the equivalent beam and the framed structure. The comparison in terms of absolute, differential settlements and angular distortion (Fig. 13) reveals a negligible role of the superstructure inertia, with variation between simplified and complete calculations ranging between − 4% and + 7% for all examined cases. This outcome was also seen with numerical calculation by Karamitros et al. (2013b), who gave a deviation lower than ± 5%. Both observations combine with the field evidence gained by Yoshida et al. (2001) in Adapazarı (Turkey) during the 1999 Kocaeli earthquake, where a limited influence of the height was noticed on the buildings damaged by liquefaction. On the contrary, a relevant role of height was seen on the buildings damaged by shaking in the non-liquefied areas of the city. This result has been explained by Karamitros et al. (2013b) considering that the seismic isolation-induced at the upper levels by the liquefied layer (see also "Appendix 1"), inhibits building oscillation and reduces inertial effects.

Prediction of absolute settlement with the Bullock et al. (2018) formula
In the last decade, several attempts have been made to produce predictive formulas to massively estimate the performance of buildings distributed over a territory. With this aim, the recent literature proposes a variety of relations to predict liquefaction-induced settlements based on the main characteristics of the phenomenon, i.e. seismic input, subsoil properties and foundation-bearing pressure (Dashti and Bray 2013;Karamitros et al. 2013;Bray and Macedo 2017;Bullock et al. 2018). Among them, the physics-informed semi-empirical probabilistic formula proposed by Bullock et al. (2018) is the last and probably the most comprehensive. The method, calibrated with the results of a 3D fully coupled numerical parametric study and validated with a database of field observations and centrifuge experiments, estimates the expected mean foundation settlement for a given set of input parameters. In this study, this formula has been applied to each analyzed condition comparing the results of numerical calculation with the value computed with the spreadsheet provided by the authors (https:// shide hdash ti. com/ geote ch-links/). The latter calculation requires inputting the parameters listed in Table 4. The cumulative absolute velocity ( CAV ) has been computed for each combination of seismic events (Table 2) and scaling factor, implementing the equation proposed by Bullock et al. (2017). Analogously, foundation contact pressure (q), width (B), depth of embedment (D f ), thickness of the non-susceptible crust (H C ), liquefaction-susceptible layer (H L ), depth from the foundation base to the centre of the liquefaction-susceptible layer (D L ) have been set equal to those given in the numerical calculation. The number of storeys (N st ) has been derived as a function of the bearing contact pressure, considering a ratio q/ N st ≈25 kPa. Finally, the mean CPT normalized cone tip resistance (q c1N ) of the liquefiable sandy layer has been computed with the equation proposed by Boulanger and Idriss (2014), as a function of the relative density (D r ): The Bullock et al. (2018) method provides a cumulated probability curve of settlements having Log-normal distribution and standard deviation σ ln equal to 0.67, these properties (2) q c1N = 0.9 ⎛ ⎜ ⎜ ⎝

Prediction of angular distortion
The analysis of the angular distortion has been carried out by looking at two classical studies on this subject for static loading, produced by Boscardin and Cording (1989) and Grant et al. (1974). The former authors proposed a vulnerability criterion applicable to brick bearing walls and small framed structures subjected to various underground events, from shallow and deep excavation to deformation induced by building construction. Their study focuses on the combination of angular distortion and horizontal strain at the foundation level, finding that the relative composition of these two kinematic components is dictated by the disturbing factor (see Fig. 15). While underground excavation tends to produce larger horizontal strains, angular distortion is predominant in the case of foundation loading. The two quantities, computed for each simulation of the present parametric study, give an alignment of dots in Fig. 15 on the line characteristic of deformation induced by the self-weight of buildings. This outcome implies that foundation movements occur with the same pattern, whether they are caused by static loading or liquefaction, being a possible explanation that subsoil deformation occurs at relatively shallow depths in both cases. Another interesting result is obtained by reporting angular distortion and maximum settlement obtained from each numerical calculation on the bi-logarithmic plot of Fig. 16, in analogy with the analysis performed by Grant et al. (1974). These authors collected observations on buildings founded on cohesive or cohesionless soils with shallow isolated or continuous footings and inferred an envelope curve for all monitoring data given by the following equation: Interestingly, this curve reported with a continuous black line in Fig. 16 represents the upper bound also for the angular distortion computed in the present study. This result highlights once more that, despite reaching different absolute values, liquefaction-induced Fig. 15 Key evidence from the parametric study: ground horizontal strain versus angular distortion overlapped to the plot of Boscardin and Cording (1989) settlements on low-rise buildings possess similar characteristics as those induced by static loading. More particularly, this line interpolates quite closely the results obtained for nil bending stiffness. For increasing EI values, the angular distortion corresponding to the same maximum settlement tends progressively to diminish. This finding is consistent with the trend seen in Fig. 12c. Figures 15 and 16, and suggests the possibility of adopting for liquefaction the same criteria used for the vulnerability assessment upon static loading. However, to apply this strategy to the extensive assessment of building over large areas (e.g. urban systems), it is necessary to predict angular distortion with fast and reliable tools. With this goal in mind, the trends shown in the examples of Figs. 6, 7, 8, 9, 10, 11 and 12 have been systematically interpreted with an analysis of variance (ANOVA, Fisher 1918). The ANOVA allows to determine the mutual influence of variables and separate those factors possessing a statistical relevance from others producing limited or random effects. Thereafter the relation between angular distortion and these variables has been sought through training an artificial neural network. This tool has been preferred to classical inference based on assigned mathematical functions thanks to its higher flexibility. Additionally, a minimum set of variables has been identified to render prediction more practical, e.g. in the form of graphical plots. After several trials, the median settlement computed with the Bullock et al. (2018) method, thickness and undrained shear strength of the crustal cap, bending stiffness of the foundation system have been identified as those compromising the lowest number of variables with an acceptable accuracy of prediction. Table 5 reports the results of the ANOVA test, where F is computed with the Fisher distribution and p measures the probability that a correlation could be random, i.e. p-value tending to zero means a greater statistical significance of the observed relation. The very low p-values computed for the assumed variables (Table 5), all lower than the suggested threshold (α = 0.05), confirm the reliability of relations. Among variables, the median settlement computed with the Bullock et al. (2018) method has the virtue of including most of the relevant factors (see Table 4) and provides a good estimate of the mean settlement (Fig. 14). This variable has been selected also considering the dependency of angular distortion on the maximum settlement seen in Fig. 16 and the relation between mean and  Grant et al. (1974) maximum settlement seen in Fig. 5. The remaining variables (H c , s u and EI) have been chosen because the one-way ANOVA test gave a significant outcome, i.e. indicated a statistical significance of the relation between the selected variables and angular distortion.
These variables have thus been assigned as input of an artificial neural network (ANN), trained with the output of numerical calculation, and used to predict angular distortion in general cases (Fig. 17a). The ANN is a mathematical model that replicates brain functioning, composed of interconnected neurons (McCulloch and Pitts 1943). The network herein assumed is a two-layer feed-forward type with sigmoid hidden neurons and linear output neurons. The information moves only in one direction from the input nodes through the hidden nodes to the output nodes, avoiding loops or recursive programming. The Levenberg-Marquardt back-propagation algorithm has then been used for training as it minimizes the sum of squares of nonlinear functions. The developed network has 10 hidden layers and the fractions of the dataset used for training, validation and testing are respectively equal to 70%, 15%, 15%.  The prediction, evaluated in terms of mean squared error MSE = 5.98 * 10 − 6 ( Fig. 17b), reveals a good performance over the whole investigated range. The statistical analysis of error shows (Fig. 17c) a symmetric distribution well approximated with a gaussian probability density function (PDF(err)) having zero mean and standard deviation s(err) = 0.0007.
The proposed tool can be downloaded at the LAGGS-geotechnics, geology and roads laboratory website with the recommendation that prediction is reliable for values of the input variables included in the investigated range, i.e. w Bullock ≤0.3 m, EI ≤ 260 MN*m, H c ≤6 m and s u ≤150 kPa. An example of ANN outcome is provided in Fig. 18, where angular distortion b is plotted as a function of median settlement w Bullock for selected values of flexural stiffness (EI equal to 0, 65 and 269 MN*m), undrained shear strength (s u equal to 50 and 100 kPa) and crust thickness (H c equal to 2, 4 and 6 m). The plots show the predominant role of the equivalent bending stiffness of the foundation system, and the less relevant but still appreciable role of the upper crust layer.

Conclusion
The objective of the present study is to derive a simple yet comprehensive method to scrutinize buildings distributed over large areas against liquefaction risk. To this aim, the attributes of liquefaction-induced settlements on low-rise buildings have been thoroughly examined with numerical Finite Difference analyses carried out with variable seismic, subsoil and structural conditions. The back-analysis of a documented case study has revealed that the pore pressure build-up and liquefaction occurring at certain depths inhibit propagation of the seismic signal to the upper levels, and thus liquefied soils act as an impedance to the shaking on the ground surface. This effect limits the role of inertia on the liquefaction assessment of buildings and justifies the herein adopted schematization of the structure with an equivalent beam. A specific calculation has shown a few percent influence of the inertia on the absolute and differential settlements computed with the simplified scheme. The parametric study has shown that, despite reaching larger values, movements induced by liquefaction present similar characters to those induced by static loading. Like for static loading, the horizontal deformation is negligible and the angular distortion b increases with the absolute settlement (mean or maximum), this relation being enveloped by the upper bound curve experimentally observed by Grant et al. (1989). This evidence has suggested to use a vulnerability criterium for liquefaction assessment from the practice of foundation engineering under static conditions, i.e. adopting the angular distortion b as a demand variable. A learning machine model (Artificial Neural Network) has thus been proposed to quantify b as a function of the fundamental characteristics of earthquake, building and subsoil. For the sake of simplicity, the number of input variables has been minimized by grouping most variables into the median settlement estimated with a semi-empirical method recently proposed by Bullock et al. (2018). The other considered input variables, selected based on their statistical relevance on b, are the flexural stiffness of the building-foundation system, the thickness and the undrained shear strength of the crust layer. The former complies with the observation of Bird et al. (2006) who distinguished absolute and differential settlement as damage factors for rigid and flexible foundations, respectively. The latter, in accordance with Karamitros et al. (2013a), who quantified the positive contribution on settlement reduction given by the non-liquefiable crust. A weblink and some design charts have been provided to compute b (the estimate error is normally distributed with a standard s(err) = 0.0007).

Appendix 1
The prototype model has been firstly validated with the back-analysis of a building located in the Municipality of Terre del Reno, whose data are extracted from the Emilia-Romagna Regional database. In particular, the MUDE, an Italian post-earthquake survey form, was consulted to derive information related to building typology and occurred damages. The building consists of a two storeys masonry building with a rectangular layout (length 13.10 m, width 11.10 m) founded on a shallow slab of poor structural characteristics. The liquefaction induced by the May 20th earthquake (M w = 6.1) caused significant differential settlements in the East-West direction, with an absolute settlement of 35 cm on the West side and 5 cm on the East side ( Fig. A1.b). The epicenter was located about 15 km far from the building, at a depth approximately equal to 10 km below the ground level (Luzi et al. 2019). The acceleration time history assigned for the back analysis (Fig. 2.c) has been computed transferring to the considered site the signal recorded at the nearest seismic station (Mirandola, from the ITACA seismic catalogue) performing the procedure suggested by Sinatra and Foti (2015). These authors propose to deconvolve the acceleration time history recorded at the station to recover the signal to the seismic bedrock. Then apply the attenuation law proposed by Bindi et al. (2011) to move from the station to the studied site, then perform a local seismic response analysis to obtain the input at the model's base. In the analysis, the deeper subsoil model has been taken from Fioravante et al. (2013), while the top subsoil stratigraphy has been reconstructed considering various CPTU tests performed in the closest area around the building (Fig. 19a).
Considering its known structural characteristics, the building has been simulated with an elastic body of given width (B = 12 m), flexural stiffness (EI = 60 MN*m) and contact pressure (q = 35 kPa). The performed analysis returned a significant rotation of the building with a final vertical displacement at the west corner of 37 cm and at the east side corner of about 4 cm, these values in close agreement with the post-earthquake surveyed deformation (Emilia-Romagna Regional database, MUDE) giving respectively 35 and 5 cm at the two corners (Fig. 19).
The results of the simulation are presented in Fig. A2 as a typical output of calculation. Fig. A2.c reports the acceleration time histories for points A (at the foundation level), B (in the middle of the liquefiable layer) and C (at the bottom of the model). During the first 36 s, the wave motion propagates upward with a slight signal modification. After this time, the pore pressures in the sandy layer start to increase, as shown by the plot of r u = Δu/σ' zo (Fig. 20e), and this effect progressively impedes the transmission from the sandy layer to the upper foundation. Acceleration is still conveyed form the clayey stratum. When r u approaches relatively high values (say 0.8), after about 38 s, the sandy soil softens and becomes unable to further propagate motion upward. This effect proceeds until liquefaction fully develops (r u approaching 1). Therefore, the increase of pore pressures and the liquefaction-induced decay of stiffness in the sandy soil, prevents wave propagation and the liquefied soil acts as an isolator for the building shaking. From this viewpoint, this effect is positive as it should reduce structural damage. On the other hand, the increase or pore pressure in the sandy soil determines (Fig. 20.) significant settlements that start to increase when r u reaches values of about 0.6 and continues with the progress of liquefaction. The maps of Fig. A4.a and b, which report the fields of pore pressure ratio R u at the end of the shaking and of maximum shear deformation throughout the event, reveal that liquefaction does not occur uniformly in the sandy layer but starts principally near singularities of the subsoil. The settlements time histories of the foundation reported in Fig. 20c show a significant differential settlement between the left and right corners of the building, with vertical displacement respectively equal to 37 and 4 cm. This effect clearly depends on the higher thickness of the liquefiable layer on the left side of the building. In conclusion, the above results show the ability of the adopted numerical model to reproduce the phenomena taking place at the building foundation upon seismic shaking and liquefaction. Despite the calculation has been simplified (a two-dimensional model has been implemented) the results are fully consistent with the observed response of the building and, particularly, with the differential settlements experienced by the foundation at its two opposite corners. Fig. 20 Results of the simulation in terms of a pore-pressure ratio field; b maximum shear field; c acceleration time histories at the bottom of the model (point C), in the middle of the liquefiable layer (point B) and at the foundation's level (point A); d settlements time histories at the left corner of the foundation (in black) and at the right corner of the foundation (in red); e pore-pressure ratio time history at the middle of the liquefiable layer.