The possible use of equivalent homogeneous subsoil models for 1D seismic response analyses in seismic microzonation studies

The possibility is here explored to use an ‘equivalent’ homogeneous configuration to simulate 1D seismic response of heterogeneous engineering-geological bodies when relatively weak seismic impedance contrasts internal to the bodies (it was assumed a shear wave velocity variation between the alternating layers equal to 150 m/s) only exist above the seismic bedrock. This equivalent configuration is obtained by considering an equivalent Vs value the harmonic average of the actual Vs values and a linear combination of G/G0 and D curves relative to the lithotechnical components present in the actual configuration. To evaluate feasibility of this approach, a wide set of numerical simulations was carried out by randomly generating subsoil layering including sequences of alternating thin layers of geotechnical units (e.g. sands and clays) each characterized by a characteristic nonlinear curve. Outcomes of these simulations are compared with those provided by considering a single homogeneous layer characterized by equivalent nonlinear curves obtained as a weighted average of the original curves. By comparing the heterogeneous and the homogeneous columns seismic response in terms of amplification factors and fundamental period, the results confirm the possibility to model a 1D column characterized by a generic lithostratigraphic succession with an equivalent one without introducing significative errors that, at least for the studied cases, do not exceed the 6%. This conclusion is substantially confirmed by extending the comparison to a real case, i.e., the 113 m-thick heterogeneous soil profile at Mirandola site (Northern Italy), presented in the last part.


Introduction
Accounting for nonlinear behavior of soil under dynamic loads is necessary to provide a consistent modelling of local seismic response during strong earthquakes (Kramer 1996). Two methods of analysis can essentially be chosen to this purpose: equivalent linear, based on a series of iterative analyses assuming a visco-elastic soil behavior, or a true nonlinear approach, this latter involving a broad range of simplified (i.e. hyperbolic family models) and advanced soil constitutive models (Hashash et al. 2010). The nonlinear curves, i.e. the normalized shear modulus (G/G 0 -γ) and the damping ratio (D-γ) curves, are directly employed in the equivalent linear strategy or used to calibrate the constitutive parameters in the true nonlinear approach. In the Seismic Microzonation (SM) studies the equivalent linear analysis is still the preferred approach because the limited data generally available on cyclic behavior as well as the large number of numerical simulations to be carried out (Pagliaroli 2018).
A crucial step in these analyses is defining variation of shear modulus and damping ratio when excitation strength increases (G/G 0 -γ and D-γ curves). These can be obtained ad hoc from laboratory cyclic or dynamic tests (e.g., resonant column, cyclic simple shear, torsional shear tests). In general, these data are provided by testing relatively small homogenous specimens assumed to be representative of the geological body of interest. However, actual soil structure is characterized by significant lithological heterogeneities (e.g. lithotypes alternating with depth) at the scale of meters or less in the lack of sharp seismic impedance contrasts. This is typically the case of flyshoid formations characteristic of many recent tectonic belts, of wide sedimentary basins affected by glaciation cycles, coastal plains in the presence of ingression/regression cycles, etc.
Therefore, in these particular cases when the lithological heterogeneities are significant, a detailed characterization representative of each lithological units, in most cases, cannot be afforded due to relatively high costs of this kind of analyses. Moreover, technical difficulties also prevent an extensive characterization of coarse grained soils (sands, gravels, pebbles, etc.). This last problem, also impose the use of modulus reduction and damping available in the technical literature and representative of data obtained worldwide where expensive studies were allowed to be performed.
The above problems are also enhanced when extensive seismic response studies are of concern and detailed geotechnical soil characterization cannot be easily performed. This is the case of site response analyses for SM, where subsoil investigations are carried out at urban scale and therefore a detailed characterization of stratigraphic profile as well as mechanical proprieties is not available in the numerical simulations. As an example, no more than two undisturbed samples for each municipality were available in the SM studies for the reconstruction after the Central Italy 2016 seismic sequence (Pergalani et al. 2020). Similar conditions can be found also in Moscatelli et al. 2012 and in other not Italian SM sites as in the case of Giallini et al. 2021;Mancini et al. 2021;Tsereteli et al. 2021. Another problem arising in these studies is the general lack of borehole data allowing a direct lithotechnical characterization of the soil structure. In most cases (in particular in the Italian practice), soil structure is inferred by geological considerations and indirect geophysical investigations (e.g. SM-WG 2015; Albarello et al. 2015;Albarello 2017;Moscatelli et al. 2020). This implies that soil characterization is performed in terms of engineering-geological categories, (TCSM 2020(TCSM , 2021Amanti et al. 2020), by grouping the different rock-units into two main categories: "Cover terrains" and "Geological bedrock" units, on the base of their lithological feature and stratigraphic position. Within each of these wide groups, a number of subgroups (Engineering Geological Units, EGU in the following) is identified by considering macroscopic features and other geological considerations (e.g. Catalano et al. 2019;Amanti et al. 2020). These EGU, are inherently heterogenous in terms geotechnical behavior and this makes difficult their use in seismic response modelling.
For SM purposes, an attractive simplification is to replace the heterogeneous subsoil model with a homogeneous equivalent one characterized by a very close seismic response at surface (in terms of response spectra and amplification factors). This operation requires the definition of an equivalent shear wave velocity, Vs ,eq , and equivalent variation curves of the shear modulus G eq /G 0 -γ and the damping ratio D eq -γ.
In this paper, a methodology was proposed to test the possibility of define homogeneous columns equivalent to the heterogeneous ones actually representative of the EGU unit of concern. For the scope, an extensive parametric study, consisting in 1D frequency-domain equivalent linear site response analyses, was carried out.
Overall, in presence of lithostratigraphic alternations like the ones studied in this work (i.e., sand-silt and sand-gravel alternations) grossly characterized in terms of relative importance of each lithotype within the stratigraphic succession, the possibility to consider in computations an 'equivalent' homogeneous counterpart defined by only considering the percentage of the soil column represented by each lithotype, seems to be feasible, at least when relatively weak impedance contrasts (i.e., ΔV s ≤150 m/s) only exist between the internal alternating lithologies, overlaying the seismic bedrock.
In the first part of the paper, some preliminary analyses are presented aiming at testing a possible strategy to define equivalent curves representative of the nonlinear behavior of a heterogeneous body. In particular, some simple configurations are considered on purpose. In the second part, outcomes of an extensive analysis are presented by supporting effectiveness of the proposed approach finally tested with the real case of Mirandola site in north of Italy.

Preliminary analysis
Consider a 1D heterogeneous column where the upper layer is constituted by sand and the lower one by silt overlaying the bedrock. Both layers have the same thickness (H 1 = H 2 = 25 m), the same shear wave velocity (Vs 1 = Vs 2 = 350 m/s), the same unit weight (γ 1 = γ 2 = 20 kN/m 3 ) while differ for the nonlinear response since different curves of the normalized shear modulus, G/G 0 , and damping ratio, D, with shear strain, γ, were assigned: for sands we considered Seed and Idriss (1970) mean curve, for silts, Vucetic and Dobry (1991) curves relative to Plasticity Index (PI) 15-30-50. The bedrock is characterized by the same unit weight and Vs of the overlying layers in order to avoid the presence of seismic impedance contrasts at the bottom of the model.
In order to define a 1D equivalent column (i.e., model characterized by the same seismic response of the heterogeneous column) representative of the corresponding geologicaltechnical unit as defined within the SM studies, it is necessary to calibrate equivalent variation curves (G eq /G 0 -γ and D eq -γ) able to provide, on average, the seismic response of the target column at surface. The input motion is a sinusoidal signal (frequency equal to 5 Hz and number of cycles equal to 10) with amplitude ranging from 0.001 to 2 g.
The procedure employed to set equivalent V s and nonlinear curves are described in the following.

3
The equivalent normalized shear modulus curve was calculated by means of the following equation where Vs ,eq is the ratio between the total column height H and the travel time, dt, of the wave propagating from the bedrock up to the surface of the heterogeneous model. The relevant dt values are obtained by considering phase shift between input (time history at the bottom) and output (time history at the free surface) signals. Repeating this procedure for all considered wave amplitudes (from 0.001 g to 2 g) allowing to explore the response at increasing shear strains, the corresponding values of Vs eq and G eq are shown in Table 1, for Vucetic and Dobry (1991) PI = 30 curves for the silt layer.
For the construction of G eq /G 0 -γ curve, the shear strain levels corresponding to the different amplitudes of input motion were obtained as harmonic mean value of the maximum shear strain profiles from equation as shown in Fig. 1a and finally multiplied for 0.65 to obtain the effective value. For the construction of D eq -γ curve, the following equation was obtained by inverting the amplification function expression available for the case of homogeneous visco-elastic layer overlaying deformable substratum and was used to calculate the equivalent damping associated to the corresponding stain level. In Eq. 3, A(f 1 ) is the amplitude of the amplification function in correspondence of the first natural resonance frequency while I is the contrast of impedance between the bedrock and the equivalent shear wave velocity.
(1) Table 1 Example of G eq /G 0 -γ and D eq -γ calibration parameters of the equivalent column assuming for the heterogeneous column Seed and Idriss (1970)  Note that impedance contrast I is higher than 1 in the nonlinear range due to the decay of shear modulus at increasing of shear strain.
In this way, equivalent decay curves G eq /G 0 -γ and D eq -γ have been obtained (Fig. 2). All parameters calculated to define G eq /G 0 -γ and D eq -γ are reported in Table 1 considering, as example, only Vucetic and Dobry (1991) PI = 30 curves for silt. The equivalent decay curves are compared with the corresponding target variation curves of the normalized shear modulus, G/G 0 , and damping ratio, D, with shear strain, γ, adopted for sand (i.e. Seed and Idriss (1970) mean curve) and for silt (i.e. Vucetic and Dobry (1991) PI = 15-30-50 Fig. 1 Outcomes of a 1D linear equivalent numerical model for the heterogeneous column including sand and silt by considering different dynamic loads. As concerns nonlinear behavior of sands the mean curves provided by Seed and Idriss (1970) were considered. The curves provided by Vucetic and Dobry (1991) for PI = 30 have been adopted curves for silt. a Maximum shear strain profiles and b amplification functions for signal amplitudes ranging from 0.001 g to 2 g Fig. 2 Comparison of equivalent the normalized shear modulus and damping ratio curves with the literature curves adopted in the corresponding heterogeneous columns curves) layers in the target column (Fig. 2). As expected, the comparison shows that the equivalent decay curves are arranged on average in the median position with respect to the literature curves adopted for the heterogeneous column, corroborating the hypothesis that, in presence of lithostratigraphic successions, even more complex than that studied in the example, that are very common in the geological-technical units defined within the SM studies, the nonlinear variation of the shear modulus and the damping ratio could be modelled with good approximation as average value weighted respect the abundance of the considered lithotype within the whole column.
A simpler procedure to compute equivalent properties has been then explored with reference to the heterogeneous columns 1 and 2 schematized in Fig. 3: both heterogeneous columns were modelled, implementing different sand-silt alternations at depth and assuming the same sand and silt percentage equal to 50%. The sand lithotype was modelled assuming Vs sand = 350 m/s, γ sand = 20kN/m 3 , Seed and Idriss (1970) mean curve, the silt lithotype was modelled assuming Vs silt = 350 m/s, γ silt = 20kN/m 3 , Vucetic and Dobry (1991) PI = 30; the seismic bedrock was modelled assuming Vs bedrock = 800 m/s, γ bedrock = 20kN/m 3 and D bedrock = 1%. The seismic response of these heterogeneous columns, subjected to four real earthquakes shown in Fig. 3 selected for increased Peak Ground Acceleration (PGA) values and with different frequency content, was compared to the response of equivalent column, subjected to the same input signals and characterized by an equivalent (Vs) eq_mean , computed by equation (equal to 350 m/s in this case), and (G/G 0 ) eq -γ and (D eq )-γ computed as a weighted average by means the equations Please note that the equivalent shear wave velocity calculated through the Eq. 4 is very close to.
Vs ,eq calculated in Eq. 1 as ratio between the total column height H and the travel time, dt, as reported in the Table 1.
In order to compare the difference between heterogeneous and equivalent models with uncertainties associated to the approach to include soil nonlinearity, linear equivalent and true nonlinear site response have been executed by means of DEEPSOIL software (Hashash et al. 2017): Fig. 4a and Fig. 4b in particular, show the comparison between nonlinear and equivalent linear analysis for the heterogeneous columns and the corresponding equivalent ones, for the target cases 1 and 2, respectively.
The comparison was proposed in terms of amplification effects. The ground motion is usually expressed in terms of peak parameters or, more frequently, integral parameters (Pagliaroli 2018), such as the integral of spectral accelerations or pseudo-velocities in a prescribed range of period. In particular, in SM studies carried out in Italy, amplification effects are generally quantified in terms of the Amplification Factor (AF) over three intervals of vibration periods in the form (Pergalani et al. 2020), where the subscripts i and o refer to the reference (input at outcropping bedrock) and soil surface (output) spectra, respectively. The AFs are represented by histograms of the average values of the four selected inputs in Fig. 4. As a whole, it can be observed a good match between the seismic responses of heterogeneous and the equivalent columns, the differences are in the same order of the uncertainties associated to the analysis method, at least for the considered signals: the equivalent linear approach slightly overestimates the AFs with respect to the nonlinear approach, particularly at low and medium periods, in any case less than 5%, while the difference between the heterogeneous and the equivalent column is contained within 6%. Furthermore, it should be noted that the maximum shear stress (γ max ) reached by means of nonlinear analyses does not exceed 2% for the higher PGA value assumed for the analyses; this result is congruent with the difference that we found between the homogeneous equivalent and the heterogeneous column and, at the same time, justifies the choice to adopt the linear equivalent analysis approach in the following parametric analyses (Kaklamanos et al. 2013).
The results obtained from this preliminary analysis confirms the possibility to model a 1D column characterized by a generic lithostratigraphic succession with an equivalent one representative of a geological-technical unit within the SM without introducing significative errors.
However, to get more general conclusions, we need to extend this study to a significant number of cases to take into account the high vertical variability typical of the subsoil conditions (Rathje et al. 2010;Kaklamanos et al. 2013;Pagliaroli et al. 2014a,b;Kim et al. 2016;Foti et al. 2019;Fabozzi et al. 2019aFabozzi et al. , b, 2021. For the scope, an extensive parametric analysis is discussed hereinafter (Sects. 3 and 4), with reference to the cases of sand-silt and sand-gravel 1D lithostratigraphic successions.
Note that the considered lithostratigraphic successions, with reference to the EGU defined within the SM standards (TCSM 2020, 2021. See Table 2), include most of the cover terrains with exception only of few units (i.e., OL, OH, CH, PT).

Parametric analysis
The parametric study consists in 1D frequency-domain equivalent linear site response analyses carried out by means of numerical simulations. Simulations concern randomly variable vertical lithostratigraphic configurations as conceptually schematized in Fig. 5. For Table 2 Engineering geological units from SM standards (TCSM 2020(TCSM , 2021 EGU from SM GW: Well sorted gravels, mixed gravels and sands GP: Not sorted gravels, mixed gravels and sands GM: Silty gravels, mixed gravels, sands and silts GC: Clayey gravels, mixed gravels, sands and clays SW: Well sorted sands, mixed sands and gravels SP: Not well sorted sands SM: Silty sands, mixed sands and silts SC: Clayey sands, mixed sands and clays OL: Organic silts, low plasticity organic silty-clays OH: Middle plasticity organic clays, organic silts MH: Inorganic silts, fine sands, diatomic silts ML: Inorganic silts, fine silty-clayey sands, low plasticity clayey, silts CL: Middle-low plasticity inorganic clays, gravel-sandy clays, silty clays CH: High plasticity inorganic clays PT: Peat and organic soils a generic lithostratigraphic succession including more than one lithotype, the column is discretized into brick elements (step 2 in Fig. 5) by assigning their maximum dimension (10 brick elements 3 m thick in the example in Fig. 5). Subsequently, all the possible bricks combinations, defined hereinafter as 'permutations', are generated (step 3 in Fig. 5). In particular, starting from the percentage of each lithotype in the lithostratigraphic succession (equal to 10% for the lithotype 1 and equal to 90% for the lithotype 2 in the example in Fig. 5), all possible not repeated permuted lithostratigraphic successions are simulated by means of the multinomial coefficient expressed as where n is the total number of bricks and K is the vector representing the bricks number for each considered lithotype.
Tables 3 and 4 synthesize the main properties assigned to the considered combinations. In both cases, 9 combinations were considered, assuming the variation of sand, silt and gravel ranging from 90 to 10%, setting the shear wave velocity, Vs i , equal to 250 m/s, 400 m/s and 550 m/s for silt, sand and gravel, respectively. Note that in the last columns of Tables 3 and 4 the number of permutations computed for each percentage combination are shown: the maximum number of permutations, equal to 252, is reached when sand, silt and gravel are equal to 50%, while the minimum number of (8) permutations, equal to 10, is reached when the relative abundance of lithotypes is 90%-10% or 10%-90%. As example, Fig. 6 shows some computed permutations for the combination 5# (sand 50%-silt 50%) including complex Vs profiles generated with multiple alternations along the depth. Regarding the nonlinear curves, the variation of the normalized shear modulus, G/ G 0 , and damping ratio, D, with shear strain, γ, has been defined according to the curves by Stokoe et al. (2004) for well-graded gravel (GW), well-graded sand (SW) and lowplasticity clay (CL) resulting from a large number of laboratory tests. As it can be observed in Fig. 7a-a', the selected curves clearly identify the different behavior of the three lithotypes, both in terms of shear modulus decay and damping increment. Other "standard" literature curves do exist: the Stokoe et al. (2004) SW and CL curves well match the well-known Seed and Idriss mean (1970) and Vucetic and Dobry IP = 15 (1991) ones, respectively; on the contrary Stokoe et al. (2004) GW curve significantly differs from Rollins et al. (1998) curves, frequently used for gravelly soils, that is more linear and very close to Seed and Idriss mean (1970). For our purposes, to assume a more pronounced difference between the cyclic behavior of the two lithotypes (sands and gravels), we preferred Stokoe et al. (2004) GW for gravelly soils.
These selected curves were adapted to all equivalent columns representative for each permutation, weighing the percentage contribution of the three considered lithotypes by means of Eqs. 5 and 6. Figure 7b-b' and Fig. 7c-c' shows the so calculated equivalent (G/G 0 ) eq -γ and (D eq )-γ curves for sand-silt and sand-gravel lithostratigraphic alternation, respectively. The seismic bedrock instead, was modelled assuming Vs bedrock = 800 m/s, γ bedrock = 20kN/m 3 and D bedrock = 1%.
Finally, regarding the input signals, three groups of 10 unscaled real accelerograms, for a total of 30 real earthquakes, recorded at outcropping seismic bedrock were selected for three different PGA levels: Class 1 (PGA ≤ 0.1 g), Class 2 (0.1 g < PGA < 0.2 g), Class 3 (0.2 g < PGA < 0.4 g), respectively. The natural accelerograms were extracted from the Italian database ITACA (http:// itaca. mi. ingv. it/), Vs (m/s) Fig. 6 Some example of permutations (5 over 252) generated into the alternation sand 50%-silt 50% considering as selection criteria: (i) the magnitude ranging from 5 to 7; (ii) the epicentral distance ranging from 10 to 80 km; (iii) seismic stations located in correspondence of site class A (based on geophysical tests) or A* (based on surface geology) according to EuroCode8 (outcropping rock conditions). Tables 5, 6, 7 show in detail the main characteristics of the selected accelerograms.
As it can be observed in Fig. 8a-c, the considered PGA values are almost uniformly distributed in the three classes, from the lower to the higher extreme of each interval.

Discussion of the results
The outcomes of the parametric analysis were processed in terms of AFs (calculated through the Eq. 7) and the fundamental periods (T 0 ) of the permutated columns compared with the same parameters associated to equivalent columns. As example, these results are shown in Fig. 9 for sand (50%)-silt (50%) and sand (50%)-gravel (50%) alternations. In Fig. 9a, c, the AFs are represented by histograms separately for the three intervals of periods (i.e., [0.1 s -0.5 s]; [0.4 s -0.8 s]; [0.7 s -1.1 s]) and the three PGA classes (from light blue to dark blue for the permuted columns and from light to dark orange for the equivalent columns). Each histogram represents the mean value over all the signals falling into the considered PGA class and limited to the permuted columns, they were also averaged with respect to the number of permutations associated with the different percentage combinations between the considered lithotypes (see permutations number in Tables 3 and 4). Analogously, in Figs. 9b, d, the fundamental period of the permuted (orange points) and the equivalent (blue points) columns are compared, once again averaging the values respect the input signals of the PGA classes and the number of permutations. Furthermore, the linear elastic fundamental period calculated as T el = 4H/Vs eq is also displayed for comparison (red line).
Due to the nonlinearity, it can be appreciated an increasing trend of the period when the PGA increases. For sand-silt alternations (Fig. 9a,b), T 0 falls inside the [0.4 s-0.8 s] range of period thus justifying the higher amplifications observed in this range. Conversely, for sand-gravel alternations (Fig. 9c, d), for PGA class 1, most cases fall inside the low interval of period while for the PGA class 2 and 3 most cases fall inside the medium and high intervals of period. Accordingly, the higher AFs are in [0.1 s-0.5 s] for PGA class 1 and in  Overall, a satisfactory agreement between the response of the permuted columns and the corresponding equivalent ones is observed, both in terms of AF and T 0. The difference between the permuted and the equivalent columns is quantified in terms of ΔAF ratio and ΔT 0ratio as expressed by the equations respectively. In particular, Tables 8 and 9 show this percentage ratio for sand-silt and sandgravel alternations respectively. In the first case (Table 8), it can be appreciated a difference that does not exceed the 6% both in terms of AF and T 0 while, in the second case (Table 9), this difference is reduced by one order of magnitude (differences lower than 0.6%) due to the higher stiffness of the whole column in presence of gravel leading to a less pronounced nonlinear behavior.
A further comparison is finally proposed in Figs. 10 and 11. In function of the PGA classes and separately for the three ranges of period, the FA and T 0 of the permuted columns (colorful points) are compared with the limit cases in which the permuted columns show a percentage ratio of sand, silt and gravel equal to 100% (x pointers). These limit cases perfectly bound the variation interval of AF an T 0 of all permutations. Note T 0,el = 4H/Vs eq Equivalent column Mean permutaƟon Fig. 9 Comparison between permuted and equivalent columns in terms of (a,c) AF and (b,d) T 0 for sand (50%)-silt (50%) and sand (50%)-gravel (50%) alternations that the results shown in Figs. 10 and 11 are averaged with respect to the number of permutations and the input signals of the considered PGA class.

Case study
An application of the above proposed methodology is here discussed for the site of Mirandola (in Emilia Romagna Region, northern Italy) located on the alluvial deposits of the Po plain; the village suffered significant damages during the 2012 Emilia seismic sequence (Tertulliani et al. 2012). Figure 12a shows the simplified stratigraphic log deducted from Garofalo et al. (2016), characterized mainly by alluvial deposits with alternating sequences of silty-clayey layers of alluvial plain and sandy horizons, for a total depth of about 110 m, overlaying the seismic bedrock consisting of marine and transitional deposits of lower-middle Pleistocene age. The stratigraphic log includes also the SM classification where four main geological units were identified with the following percentages: MH = 18%; CL = 21%; SP = 19%; SW = 42%. The adopted Vs profile shown in Fig. 12b, is consistent with the corresponding profiles obtained from down-hole and cross-hole surveys performed, within the InterPACIFIC project (Garofalo et al. 2016) in areas where the Marnoso-Arenacea formation outcrops. According with the implemented methodology, Stokoe et al. (2004) SW and CL curves were adopted to model the nonlinearity of SP-SW and MH-CL geological units, in the percentage of 61% and 39%, respectively. The site response of the so defined column was compared with the corresponding equivalent column characterized by Vs, eq = 304 m/s and equivalent decay curves weighted respect the considered alternation. Seven input signals recorded on rock outcrop were extracted from the European Strong Motion database (ESMd, https:// esm. mi. ingv. it), and scaled up to the reference value of the peak ground acceleration (PGA ref ). In particular, according to the  probabilistic seismic hazard approach (PSHA), adopted by the Italian Building code (NTC 2018), the median hazard curve (50° percentile) provides a value of a g about equal to 0.13 g, in correspondence of a return period of 475y. The spectrum-compatibility of the selected signals was achieved by adopting In-Spector REV007 software (Acunzo et al. 2014) in the period range 0.1-1.1 s. The procedure required a manual pre-selection of the signals based on the disaggregation histogram of the site; in particular, the selection was set in the ranges of 4 < Mw < 7 and 0 < Ep.Dist < 60 km, encompassing the whole de-aggregation histogram. In Table 10 are synthesized the main characteristics of the selected accelerograms which average spectrum resulted compatible with the corresponding uniform target spectrum 50° percentile as shown in Fig. 12c. By comparing the seismic response of Mirandola site between the real (or target) and the equivalent column, it can be observed a good agreement in terms of Fourier Amplitude Spectrum (FAS) and AF as shown in Fig. 13a, b. Note that the spectra and the histograms in Fig. 13 are the mean value calculated respect the seven reference input signals.
∆AF ratio,mean was also calculated for the three intervals of period: ∆AF [0.1 s-0.5 s] = 6.5%; ∆AF [0.4 s-0.8 s] = 7%; ∆AF [0.7 s-1.1 s] = 6%. Please also note that the equivalent model provides conservative estimate of the amplification factors as also observed in most sand-silt alternations in parametric analysis (see Table 8). In terms of period instead, it was calculated a  value of T 0,mean for the real column equal to 0,63 s and for the equivalent column equal to 0.6 s that returned a difference between the models in terms of ∆T 0,mean about equal to 5%. The resulting comparison shows that, although for the case of Mirandola it was considered a deeper column (113 m deep) respect the column depth assumed for the parametric analyses (30 m deep), the order of magnitude of the approximation found by modelling the Mirandola real column with the equivalent one is almost the same of the parametric study for the case of sand-silt alternation (Sect. 4). Overall, this finding from one hand strengthens the proposal to model the lithostratigraphic alternation with equivalent homogeneous model and, from another hand, demonstrates the possibility to extend this methodology also to column rather deep as for the case of Mirandola site.

Conclusions
In the present work, the possibility explored to characterize nonlinear soil behavior of heterogeneous deposits by considering 'equivalent curves' deduced by combining curves relative to the lithotechnical components of the soil deposits. This possibility could be of main importance when fast-and-dirty seismic response evaluations must be extensively determined for heterogeneous soil columns (e.g., lithotypes alternating with depth), which is the case of seismic microzonation studies generally finalized to seismic urban planning, in the lack of any detailed geotechnical characterization of each soil component.
The proposal consists in replacing the heterogeneous subsoil model with a homogeneous equivalent one, representative of the Engineering Geological Units (EGU) encountered in the deposit, characterized by the same seismic response at surface (in terms of response spectra and amplification factors). This operation requires the definition of an equivalent shear wave velocity, Vs ,eq , and equivalent nonlinear curves of the normalized shear modulus G eq /G 0 -γ and damping ratio D eq -γ. The equivalent soil column is characterized in terms of curves corresponding to the weighted average of curves representative of the single components, which is function of its percentual presence.
The outcomes of a wide set of numerical 1D site response analyses, carried for the scope by randomly generating subsoil layering. These simulations, including sequences of alternating thin layers of geotechnical units that cover most of EGU, demonstrate that equivalent homogenous columns and their heterogeneous counterparts provide a seismic response (in terms of amplification and main resonance period) differing not more than 6% (see Tables 8 and 9). This discrepancy is of the same order of magnitude of that existing when different numerical modelling procedures are applied.
The highest mismatch between the target and the equivalent column was found in sandsilt alternations: this difference is reduced by one order of magnitude (differences lower than 0.6%) in the case of sand-gravel alternation due to the higher stiffness of the whole column in presence of gravel leading to a less pronounced nonlinear behavior.
The parametric study discussed above (Sect. 4) presents two main limitations in terms of maximum investigated column depth (the maximum depth was set equal to 30 m in the parametric analysis) and in terms of minimum layer thickness (the minimum layer thickness was set equal to 3 m in the parametric analysis). In order to overcome the limitation related to the maximum investigated column depth, the case study of Mirandola site characterized by a 113 m-thick soil deposit is presented in §5. The results show that, also for columns deeper that 30 m, the mismatch found by substituting the lithostratigraphic succession with the equivalent column is close to the one found in the parametric study proposed in this work, at least for the case of sand-silt alternation (Sect. 4). This could confirm the possibility to extend this methodology to a wide range of cases.
It remains instead the limitation in terms of minimum layer thickness, although we believe that this hypothesis could not affect very much the forecast of the methodology, also in consideration of the fact that generally, within the Italian SM studies, the degree of detail reached in the description of the stratigraphic successions does not exceed 3 m-thick, as considered in the initial assumption.
Finally, as already mentioned in Sect. 1, the applicability of the proposed methodology is limited to lithostratigraphic alternations exhibiting internal contrast of impedance comparable with ones considered in the present work ( ≤ 150 m/s). Thus, it is not applicable in presence of higher contrast of impedance as, for example, in the case of marine and fluvial terraces, alluvial and coastal plains, intermountain basins, cliffs and lava rock layers where the contrast of impedance plays a not negligible rule as observed in Panzera et al. 2011;Panzera et al. 2013;Imposta et al. 2018;Fabozzi et al. 2021.