One-Dimensional Nonlinear Seismic Response Analysis Using Strength-Controlled Constitutive Models : The Case of the Leaning Tower of Pisa ’ s Subsoil

The Leaning Tower of Pisa was built between 1173 and 1360 and began to lean at the beginning of its construction. Extensive investigations to reveal the causes of the tilting only began in the early 20th century. Although few earthquakes have been recorded, there is a renewed interest in the seismic behavior of the tower triggered by the availability of new data and technologies. This paper highlights the influence of using new strength-controlled constitutive models in case of 1D nonlinear response analysis. This is an aspect that has been poorly investigated. Most of the computer codes currently available for nonlinear seismic response analysis (SRA) of soil use constitutive models able to capture small-strain behavior, but the large-strain shear strength is left uncontrolled. This can significantly affect the assessment of a 1-D response analysis and the Leaning Tower’s subsoil can be useful for this study as it represents a well-documented and well-characterized site. After a geological and geotechnical description of the subsoil profile and a synthesis of available data, the seismic input is defined. One-dimensional SRAs were carried out by means of a computer code which considers an equivalent-linear soil modelling and two codes which assume nonlinear soil response and permit to use strength-controlled constitutive models. All the parameters were calibrated on the basis of the same soil data, therefore allowing for a direct comparison of the results.


Introduction
The Leaning Tower of Pisa was built between 1173 and 1360, and many attempts were made during the construction to correct the undesired tilt.The tower is in the form of a hollow cylinder surrounded by six loggias with columns and vaults merging from the base cylinder and surmounted by a belfry.The structure is subdivided into eight levels called 'orders'.The outer surfaces are made with marble, while the inner ones with various masonry materials.The annulus between the outer and inner surfaces is filled with rubble and mortar.A spiral staircase winds up within the annulus up to the 6th order, while two shorter spiral staircases lead to the floor and top of the belfry.The staircase forms a large opening on the south side just above the level of the first cornice, where the thickness of the masonry suddenly decreases.The high stress within this region was a major cause of concern since it could give rise to an abrupt brittle failure of the masonry.
Extensive instrumental measurements and investigations started only in the 20th century, as illustrated in [1].The plane of maximum inclination is approximately coincident with the north-south plane.From 1993 to 2001 the Committee for the Safeguard of Leaning Tower of Pisa carried out several interventions aimed at geotechnical stabilization [2] and structural strengthening [3].
A comprehensive description of the history of the construction of the Tower and interventions is described in [2].A dynamic monitoring system was installed by the Committee for the Safeguard of the Leaning Tower of Pisa at the end of its activities in 2001.During the last 20 years, several low-intensity earthquakes have been recorded.Although few events and very slight damage related to earthquakes has been observed from the construction of the monument, there is a considerable interest in the seismic behavior of the tower [4][5][6][7] triggered by the availability of new data and technologies.Within the chain starting from seismogenic zones and up to the tower, the dynamic response of its subsoil holds a crucial importance, particularly in terms of the frequency content of seismic input.This paper presents some results showing the influence of using new strength-controlled constitutive models in case of 1D nonlinear response analysis.This is an aspect poorly investigated up to now.Most of the codes currently available for nonlinear seismic response analysis use constitutive models able to capture small-strain behavior, but the large-strain shear strength is not controlled.The development of hyperbolic models (for example the Modified Kondner-Zelasko [8]), in the last decades permitted to well capture the backbone stress-strain curve and the unloading-reloading response of the soil.These models are fit to the experimental or reference shear modulus reduction curves and damping ratio curves but are unable to give realistic results in case of medium-large strain levels, resulting in soil layers strength underestimation or overestimation.This fact can significantly affect the assessment of a 1D site response analysis.

Geological Description of the Study Area
The current configuration of the Pisa plain comes from the dynamic interaction between erosive and depositional processes developed, starting from the end of the Last Glacial Maximum expansion, first in coastal and transitional marine environments.Then, in the Middle Ages, it evolved mainly in continental fluvial-lake environments.
The proximity to the coast provided the environmental conditions, since prehistoric times, to the development of the anthropic settlement in the area on which Pisa will rise and in the entire coastal area.Figure 1 shows a geological map of the valley of the Arno river, from Florence to the Tyrrhenian (Ligurian) Sea.The mountains positioned to the east of Pisa consists of formations belonging to the Paleozoic and the Tertiary geological periods and their structures reveal intense tectonic deformations.
The main structure is oriented towards NNW-SSE, with faults bordering the western part of the mountains and located 5 km northeast of Pisa.
More recent deposits, identified as quaternary continental in the geological map, developed in marine and in fluvial-lacustrine environments and deposited mainly in an estuary environment.
The valley of the Arno river crosses a mountainous territory starting from its source, in the Tuscan-Emilian Apennines, up to Florence (Figure 1).During the Last Glacial Maximum, in the period between 23,000 and 15,000 years ago, the level of the sea lowered at least 100 m and the river dug a deep valley west of Florence.This valley is now filled with alluvial materials.
Pisa rises on these alluvial deposits, at an elevation of 3-4 m above the present mean sea level.The thickness of alluvial sediments near the city is at least 300 m.In the Piazza dei Miracoli area, at the depth of about 40 m below the ground surface level can be found marine sands deposited during the Flandrian transgression.
deposit, which was formed in a period of rapid eustatic lifting of the sea level.Figure 2 shows the morphological map of the valley of the Arno river.

Geotechnical Model of the Leaning Tower of Pisa Subsoil
The soil profile under the tower is shown in Figure 3.It consists of three distinct horizons.Horizon A is about 10-m thick and primarily consists of estuarine deposits, laid down under tidal conditions with a rather erratic sequence of sandy and clayey silt layers.Typically of estuarine deposits, there are significant variations over short horizontal distances.At the bottom of horizon A there is a 2-m thick medium dense fine sand layer (i.e., upper sand).Horizon B consists primarily of marine clay, which extends down to a depth of about 40 m.It is subdivided into four distinct layers.The upper layer is soft sensitive clay, locally known as Pancone.It is underlain by an intermediate layer of stiffer clay, which in turn overlies a sand layer (the intermediate sand).The bottom layer of horizon B is a normally consolidated clay known as the lower clay.Horizon B is very uniform laterally near the tower.
Horizon C is composed of dense sand (the lower sand), which extends down to a considerable depth.Based on sample descriptions and piezocone tests, the materials to the south of the tower Above these dense sands layers, partly rearranged due to wind action, a complex mainly composed by clayey soils can be found.The clay layer close to the surface level is a clayey marine deposit, which was formed in a period of rapid eustatic lifting of the sea level.Figure 2 shows the morphological map of the valley of the Arno river.
Geosciences 2018, 8, x FOR PEER REVIEW 3 of 20 deposit, which was formed in a period of rapid eustatic lifting of the sea level.Figure 2 shows the morphological map of the valley of the Arno river.

Geotechnical Model of the Leaning Tower of Pisa Subsoil
The soil profile under the tower is shown in Figure 3.It consists of three distinct horizons.Horizon A is about 10-m thick and primarily consists of estuarine deposits, laid down under tidal conditions with a rather erratic sequence of sandy and clayey silt layers.Typically of estuarine deposits, there are significant variations over short horizontal distances.At the bottom of horizon A there is a 2-m thick medium dense fine sand layer (i.e., upper sand).Horizon B consists primarily of marine clay, which extends down to a depth of about 40 m.It is subdivided into four distinct layers.The upper layer is soft sensitive clay, locally known as Pancone.It is underlain by an intermediate layer of stiffer clay, which in turn overlies a sand layer (the intermediate sand).The bottom layer of horizon B is a normally consolidated clay known as the lower clay.Horizon B is very uniform laterally near the tower.
Horizon C is composed of dense sand (the lower sand), which extends down to a considerable depth.Based on sample descriptions and piezocone tests, the materials to the south of the tower

Geotechnical Model of the Leaning Tower of Pisa Subsoil
The soil profile under the tower is shown in Figure 3.It consists of three distinct horizons.Horizon A is about 10-m thick and primarily consists of estuarine deposits, laid down under tidal conditions with a rather erratic sequence of sandy and clayey silt layers.Typically of estuarine deposits, there are significant variations over short horizontal distances.At the bottom of horizon A there is a 2-m thick medium dense fine sand layer (i.e., upper sand).appear to be more silty and clayey than to the north, and the upper sand layer is locally thinner.This is believed to be the cause of the southward tilt of the tower [11].The water table in horizon A is found at a depth between 1 and 2 m below the ground surface; the latter has an average elevation of 3 m above mean sea level.Pumping from the lower sand has resulted in downward seepage from horizon A with a pore pressure distribution with depth which is slightly below hydrostatic (Figure 4).
Many borings beneath and around the tower show that the surface of the Pancone clay stretches out beneath the tower, from which it can be deduced that the average settlement of the monument is more than 3 m (Figure 3).
A comprehensive description of subsoil investigations carried out since the beginning of the 20th century can be found in [11], whereas the dynamic characterization of subsoil is the focus of the present paper.As reported in [12], several tests were performed from the early 1990s to determine the shear wave velocity profile Vs, namely Down Hole (DH) test [13], Cross Hole (CH) test, DH and CH tests [11]   Horizon C is composed of dense sand (the lower sand), which extends down to a considerable depth.Based on sample descriptions and piezocone tests, the materials to the south of the tower appear to be more silty and clayey than to the north, and the upper sand layer is locally thinner.This is believed to be the cause of the southward tilt of the tower [11].
The water table in horizon A is found at a depth between 1 and 2 m below the ground surface; the latter has an average elevation of 3 m above mean sea level.Pumping from the lower sand has resulted in downward seepage from horizon A with a pore pressure distribution with depth which is slightly below hydrostatic (Figure 4).appear to be more silty and clayey than to the north, and the upper sand layer is locally thinner.This is believed to be the cause of the southward tilt of the tower [11].The water table in horizon A is found at a depth between 1 and 2 m below the ground surface; the latter has an average elevation of 3 m above mean sea level.Pumping from the lower sand has resulted in downward seepage from horizon A with a pore pressure distribution with depth which is slightly below hydrostatic (Figure 4).
Many borings beneath and around the tower show that the surface of the Pancone clay stretches out beneath the tower, from which it can be deduced that the average settlement of the monument is more than 3 m (Figure 3).
A comprehensive description of subsoil investigations carried out since the beginning of the 20th century can be found in [11], whereas the dynamic characterization of subsoil is the focus of the present paper.As reported in [12], several tests were performed from the early 1990s to determine the shear wave velocity profile Vs, namely Down Hole (DH) test [13], Cross Hole (CH) test, DH and CH tests [11]   Many borings beneath and around the tower show that the surface of the Pancone clay stretches out beneath the tower, from which it can be deduced that the average settlement of the monument is more than 3 m (Figure 3).
A comprehensive description of subsoil investigations carried out since the beginning of the 20th century can be found in [11], whereas the dynamic characterization of subsoil is the focus of the present paper.
As reported in [12], several tests were performed from the early 1990s to determine the shear wave velocity profile V s , namely Down Hole (DH) test [13], Cross Hole (CH) test, DH and CH tests [11] and an SDMT test (Seismic Dilatometer Test).The results obtained from these tests show a satisfactory agreement of the V s values.All the tests reached a maximum depth of 40 m except for the CH test, which reached a depth of 65 m.None of the tests was successful in identifying the seismic bedrock.
The bedrock (i.e., V s ≥ 800 m/s) is presumed to be very deep and some attempts have been made to increase the knowledge of the subsoil below 65 m depth.An array 2D test has been exploited to obtain the shear wave velocity profile down to more than one hundred meters from ground surface [6]. Figure 5 shows the shear velocity profile down to a depth of 100 m.The array 2D test revealed the presence of a stiffer layer (i.e., V s ≈ 500 m/s) at a depth of about 100 m.Another methodology to investigate about the seismic bedrock depth by using the available H/V ratios is described in [14,15], however further studies and investigations have been planned in order to improve the knowledge of subsoil below 100 m depth.
Geosciences 2018, 8, x FOR PEER REVIEW 5 of 20 The bedrock (i.e., Vs ≥ 800 m/s) is presumed to be very deep and some attempts have been made to increase the knowledge of the subsoil below 65 m depth.An array 2D test has been exploited to obtain the shear wave velocity profile down to more than one hundred meters from ground surface [6]. Figure 5 shows the shear velocity profile down to a depth of 100 m.The array 2D test revealed the presence of a stiffer layer (i.e., Vs ≈ 500 m/s) at a depth of about 100 m.Another methodology to investigate about the seismic bedrock depth by using the available H/V ratios is described in [14,15], however further studies and investigations have been planned in order to improve the knowledge of subsoil below 100 m depth.
Many samples have been retrieved from the ground around the tower.Several of these samples have been subjected to laboratory tests to obtain the dynamic characterization of the subsoil.Further information can be found in [11].

Seismic Input
Defining the seismic input is a key step in forecasting the seismic action on a structure.Especially when the bedrock depth is defined, the assessment of the seismic hazard is usually based on well-established studies [16].
Ancient monumental masonry buildings are particularly vulnerable to vibrations and specifically to seismic actions due to many factors: in the case of the Tower of Pisa, the construction process lasted two centuries, its characteristic inclination and the several changes in the foundation system occurred during its life contributed to increase the uncertainty in the dynamic behavior.Therefore, it is necessary to define the seismic input in terms of natural accelerograms, which can be scaled in order to match a target response spectrum.In this context, the use of the design response spectrum defined by the Code (e.g., EC8) is not recommended [17], as it has no direct relationship with a specific earthquake characterized by magnitude M and source-to-site distance R, necessary to select sets of natural accelerograms.To this aim, it is appropriate, as reported in recent studies applied to archeological sites [18] to use a hybrid approach which can be briefly summarized as follows: 1. Definition of a Uniform Hazard Spectrum (UHS) obtained by means of a Probabilistic Seismic Hazard Assessment [19] using one or more Ground Motion Predictive Equation (GMPE); 2. Disaggregation of the seismic hazard to obtain the most likely combinations of M and R for a given return period (RP) [20]; 3. Definition of a Scenario Earthquake selected from the seismic catalogue, compatible with the previously defined UHS, and evaluated with the same GMPEs of point 1.Many samples have been retrieved from the ground around the tower.Several of these samples have been subjected to laboratory tests to obtain the dynamic characterization of the subsoil.Further information can be found in [11].

Seismic Input
Defining the seismic input is a key step in forecasting the seismic action on a structure.Especially when the bedrock depth is defined, the assessment of the seismic hazard is usually based on well-established studies [16].
Ancient monumental masonry buildings are particularly vulnerable to vibrations and specifically to seismic actions due to many factors: in the case of the Tower of Pisa, the construction process lasted two centuries, its characteristic inclination and the several changes in the foundation system occurred during its life contributed to increase the uncertainty in the dynamic behavior.Therefore, it is necessary to define the seismic input in terms of natural accelerograms, which can be scaled in order to match a target response spectrum.In this context, the use of the design response spectrum defined by the Code (e.g., EC8) is not recommended [17], as it has no direct relationship with a specific earthquake characterized by magnitude M and source-to-site distance R, necessary to select sets of natural accelerograms.To this aim, it is appropriate, as reported in recent studies applied to archeological sites [18] to use a hybrid approach which can be briefly summarized as follows: 1.
Definition of a Uniform Hazard Spectrum (UHS) obtained by means of a Probabilistic Seismic Hazard Assessment [19] using one or more Ground Motion Predictive Equation (GMPE);

2.
Disaggregation of the seismic hazard to obtain the most likely combinations of M and R for a given return period (RP) [20]; 3.
Definition of a Scenario Earthquake selected from the seismic catalogue, compatible with the previously defined UHS, and evaluated with the same GMPEs of point 1.
On the basis of past studies [6], the target response spectrum was defined for two return periods (RPs) of the seismic action, namely RP = 130 years (moderate earthquake) and RP = 500 years (severe earthquake).First, the disaggregation of the seismic hazard was performed, thus providing a scenario earthquake for each return period.Based on the Italian seismic catalogue CPTI15 [21], two controlling earthquakes were identified [12]: a M W 5.15 seismic event with R (epicentral distance) = 19 km was selected for an RP of 130 years (Livorno 1742) whereas a M W 5.71 earthquake with R = 21 km (Orciano Pisano 1846) was considered for an RP of 500 years.
The GMPE proposed by Akkar and Bommer [22] was then used to evaluate the target response spectra for a Ground Type B (according to the Eurocode 8 [23] and the Italian Building Code NTC-08 [24]), including the subsoil term of the equation depending on the V s,30 value of the site, which is an average shear wave velocity of the first 30 m. Seven accelerograms were selected for each return period from the European Strong Motion Database [25], considering 5 < M < 5.5 for RP = 130 years and 5.3 < M < 6.2 for RP = 500 years.The horizontal components of the selected accelerograms were scaled so that a good match was obtained between the average spectrum of the seven accelerograms of each set and the reference target spectrum in the range of natural periods 0.31-1 s.This was done to take into account the natural periods of the first two bending modes of the tower and that of the third mode (about 0.3 s), thus obtaining a scale factor (SF) for each record.Figures 6 and 7 depict the time histories of the scaled accelerograms for an RP of 130 and 500 years, respectively, as well as the scale factor used to obtain the final accelerograms.
On the basis of past studies [6], the target response spectrum was defined for two return periods (RPs) of the seismic action, namely RP = 130 years (moderate earthquake) and RP = 500 years (severe earthquake).First, the disaggregation of the seismic hazard was performed, thus providing a scenario earthquake for each return period.Based on the Italian seismic catalogue CPTI15 [21], two controlling earthquakes were identified [12]: a MW 5.15 seismic event with R (epicentral distance) = 19 km was selected for an RP of 130 years (Livorno 1742) whereas a MW 5.71 earthquake with R = 21 km (Orciano Pisano 1846) was considered for an RP of 500 years.
The GMPE proposed by Akkar and Bommer [22] was then used to evaluate the target response spectra for a Ground Type B (according to the Eurocode 8 [23] and the Italian Building Code NTC-08 [24]), including the subsoil term of the equation depending on the Vs,30 value of the site, which is an average shear wave velocity of the first 30 m. Seven accelerograms were selected for each return period from the European Strong Motion Database [25], considering 5 < M < 5.5 for RP = 130 years and 5.3 < M < 6.2 for RP = 500 years.The horizontal components of the selected accelerograms were scaled so that a good match was obtained between the average spectrum of the seven accelerograms of each set and the reference target spectrum in the range of natural periods 0.31-1 s.This was done to take into account the natural periods of the first two bending modes of the tower and that of the third mode (about 0.3 s), thus obtaining a scale factor (SF) for each record.Figures 6 and 7 depict the time histories of the scaled accelerograms for an RP of 130 and 500 years, respectively, as well as the scale factor used to obtain the final accelerograms.

Site Response Analysis
The subsoil model adopted for site response analyses is reported in Table 1.The thickness and unit weight for each lithotype were assumed on the basis of the data reported in [11].The assumed shear wave velocity profile is also based on the outcomes of the 2D seismic array test [6].However, this profile satisfactorily matches the Vs values measured by other geophysical tests in the upper part.According to the Vs profile, the seismic bedrock (Vs > 800 m/s) is not localized in the explored depth range because a Vs equal to 500 m/s is reached in the lower strata.For this reason, considering the uncertainties in the location of the seismic bedrock (Vs > 800 m/s), the accelerograms selected for the site response analysis were recorded at Ground Type B, according to Eurocode 8 definition.

Site Response Analysis
The subsoil model adopted for site response analyses is reported in Table 1.The thickness and unit weight for each lithotype were assumed on the basis of the data reported in [11].The assumed shear wave velocity profile is also based on the outcomes of the 2D seismic array test [6].However, this profile satisfactorily matches the Vs values measured by other geophysical tests in the upper part.According to the Vs profile, the seismic bedrock (Vs > 800 m/s) is not localized in the explored depth range because a Vs equal to 500 m/s is reached in the lower strata.For this reason, considering the uncertainties in the location of the seismic bedrock (Vs > 800 m/s), the accelerograms selected for the site response analysis were recorded at Ground Type B, according to Eurocode 8 definition.

Site Response Analysis
The subsoil model adopted for site response analyses is reported in Table 1.The thickness and unit weight for each lithotype were assumed on the basis of the data reported in [11].The assumed shear wave velocity profile is also based on the outcomes of the 2D seismic array test [6].However, this profile satisfactorily matches the V s values measured by other geophysical tests in the upper part.According to the V s profile, the seismic bedrock (V s > 800 m/s) is not localized in the explored depth range because a V s equal to 500 m/s is reached in the lower strata.For this reason, considering the uncertainties in the location of the seismic bedrock (V s > 800 m/s), the accelerograms selected for the site response analysis were recorded at Ground Type B, according to Eurocode 8 definition.Regarding the nonlinear properties, most of lithotypes were characterized with resonant column (RC) tests.The A1 lithotype, was characterized within this study through a recently performed resonant column (RC) test.The corresponding results are reported in Figure 9 in terms of normalized secant shear modulus (Figure 9a) and damping ratio (Figure 9b) variation as a function of the shear strain amplitude (G/G 0 -γ and D-γ curves).
Geosciences 2018, 8, x FOR PEER REVIEW 8 of 20 Regarding the nonlinear properties, most of lithotypes were characterized with resonant column (RC) tests.The A1 lithotype, was characterized within this study through a recently performed resonant column (RC) test.The corresponding results are reported in Figure 9 in terms of normalized secant shear modulus (Figure 9a) and damping ratio (Figure 9b) variation as a function of the shear strain amplitude (G/G0-γ and D-γ curves).Because of the lack of experimental data, G/G0-γ and D-γ curves that have been reported in the literature for similar soils were employed for the remaining lithotypes: MG, A2, B6, C1, C2 and C3 [26,27] (Table 1).Because of the lack of experimental data, G/G 0 -γ and D-γ curves that have been reported in the literature for similar soils were employed for the remaining lithotypes: MG, A2, B6, C1, C2 and C3 [26,27] (Table 1).These literature curves have been obtained by a huge database of samples and can quantitatively capture the dynamic behavior of the above lithotypes because are based on the grain size distribution, the plasticity index and the soil stress state.Considering the unavailability of the experimental curves, those proposed by Rollins [26] and by Darendeli [27] are a good compromise for the dynamic characterization of the MG and A2, B6, C1, C2, C3, respectively.
In order to compare the results obtained with different 1-D SRA models, analyses were carried out by means of three different computer codes: equivalent-linear code STRATA [28], nonlinear code DEEPSOIL [29][30][31] and nonlinear code ONDA [32,33].All parameters for the calibration of the models were deduced from the same group of in situ and laboratory tests.
The analyses were carried out considering a one-dimensional geometry.This choice is justified by the fact that all the main layers are horizontally stratified, and it can be reasonable to assume a geotechnical subsoil model indefinitely extended horizontally.
The seismic wave field was simplified with SH-waves vertically propagating from the bedrock to the soil deposit surface.The input motion was considered as an outcrop motion and was deconvoluted by the computer codes used in this work down to the depth of the base layer.

Seismic Response Analyses with Equivalent-Linear Models
The equivalent-linear (EL) analyses with STRATA were carried out in the frequency domain, and the amplification function was computed.The variation of the shear moduli and of the damping ratios were dependent only on the shear strain level.This means that the material damping has only a frequency independent component.The effective shear strain (γ eff ) for a given soil layer and for a given stress time history was evaluated following the suggestions in [34] (Equations ( 1) and ( 2)).
where R is the effective strain ratio, γ max is the maximum shear strain computed in a layer and M is the magnitude related with the acceleration time history considered.This quantity is representative of the time history duration (i.e., number of cycles).
The effective shear strain is needed in order to identify, by an iterative procedure based on G-γ and D-γ curves of each stratum, a strain level which is representative of a given stress time history and therefore to select appropriate values of the shear modulus G and damping ratio D.

Seismic Response Analyses with Nonlinear Models
Nonlinear analyses with ONDA and DEEPSOIL are performed in the time domain with a step-by-step integration of the equations of motion of a multi-degree of freedom lumped parameters model.For that purpose, in ONDA [32] the Wilson-θ time integration method was adopted [35], whereas in DEEPSOIL [29] the dynamic equilibrium equation is solved numerically at each time step using the Newmark-β method [36].
In ONDA, the constitutive law (τ-γ) of the initial loading stress-strain curve (backbone or skeleton curve) is represented by the Ramberg-Osgood model [37].This model consists of a not invertible relationship that depends on four parameters: α and R represent the position and the curvature of the skeleton curve, respectively, while τ max and G 0 represent the soil shear strength and the shear modulus at small strain levels, respectively.ONDA uses a strength-controlled constitutive model.The backbone curve can be written in this form: where x and y are the normalized shear strain and shear stress, respectively: x = γ γ re f (4) y = τ τ max (6) α and R were obtained by a linear interpolation of data collected from Resonant Column tests [11] for lithotypes A1, B1, B2, B3, B4, B5, B7, B8, B9 and B10, while for lithotypes MG, A2, B6, C1, C2 and C3 by a linear interpolation of data based on the modulus reduction curves reported in Table 1.For that purpose, the Ramberg-Osgood equation was rewritten as follows: Considering Equation ( 7) and performing a linear regression of the experimental and literature curves data, α and R were obtained.
In Equation ( 7), the soil shear strength (τ max ) was computed according to Equation ( 8), as proposed by Hardin and Drnevich [38], nevertheless other relationships can be used.
In Equation ( 8) σ' v0 is the vertical effective geostatic stress, K 0 is the at-rest earth pressure coefficient, ϕ is the angle of friction, and c' is the effective cohesion.The values of K 0 , ϕ and c' are reported in Table 1 for each lithotype.
DEEPSOIL can handle mainly two nonlinear constitutive laws (τ-γ) for the initial loading stress-strain curve (backbone or skeleton curve).The most-commonly used stress-strain relationship is the MKZ (Modified Kondner-Zelasko) model [8].This is a modified hyperbolic model that defines the monotonic shear stress-shear strain law coupled with unloading-reloading behavior and can properly capture the nonlinear soil response in case of small-strain values.Equations ( 9) and ( 10) are the basic backbone and the unloading-reloading relationships, respectively, for the MKZ model [8].
where γ r is a reference shear strain, β and s are dimensionless factors, τ rev and γ rev are the shear stress and shear strain at the reversal loading point, respectively.The dimensionless factors are computed by fitting the experimental data.Additional details about the available modulus reduction and damping curve fitting procedures are reported in [29].
Nevertheless, the MKZ model was found to be unable to correctly characterize the nonlinear soil behavior in case of medium-large shear strains because the shear strength at large shear strains are generally left uncontrolled.
In 2015, the General Quadratic/Hyperbolic (GQ/H) Strength-Controlled constitutive model [30,31] has been developed and implemented in DEEPSOIL code.This model can consider both the soil shear strength at failure and the small-strain stiffness nonlinearity.The GQ/H model can overcome the above limitations in the MKZ model, related to cases in which the shear strength of the soil is underestimated or overestimated leading to inaccurate site response analysis results.Further details can be found in [30,31].Compared to the MKZ model, in the GQ/H model both the initial shear stiffness (G 0 ) and the maximum shear strength (τ max ) must be defined and inputted, and the backbone curve is represented by the Equation (11).
θ τ is a curve-fitting parameter defined in the Equation (12).
In the nonlinear codes (ONDA and DEEPSOIL) used in this work the cyclic behavior of unloading-reloading (hysteretic curves) can be modelled using either the second Masing criterion [39] or the modified second Masing criterion [32,40].
The second Masing criterion states that the unloading-reloading branches of the stress-strain curve have the same shape as the skeleton curve but with a scale factor (n) equal to 2. This function can be expressed with Equations ( 10) and ( 13) for the code DEEPSOIL and the code ONDA, respectively.
where x c and y c are the normalized amplitudes of strain and stress at the stress reversal point, respectively.The hysteretic damping computed using the unloading-reloading stress-strain loops according to the Masing rules may lead to an overestimation of damping at large strains.
The modified second Masing criterion [40] assumes that the scale factor can be different from 2. In ONDA to simulate cyclic hardening behavior, the scale factor n should be higher than 2, whereas cyclic softening or material degradation could be modelled by decreasing the values of n (even n < 2).The analyses are conducted in terms of total stresses.
The n values should be determined experimentally.If experimental data are lacking, ONDA assumes that n is dependent on the over-consolidation ratio [32], and considers n equal to 5 and 3.5 in NC and OC soil conditions, respectively.When the number of cycles is increased, n decreases.
In DEEPSOIL, the modified second Masing criterion is formulated in a different way, as described in [27,31].We thus compared the site response analysis results obtained with these two codes following the standard second Masing rule.
In addition to the hysteretic damping due to the constitutive models adopted, ONDA and DEEPSOIL consider the viscous damping at small strain levels with a damping matrix expressed as a linear combination of the mass and the stiffness matrixes (i.e., Rayleigh damping formulation).Consequently, even at small strain levels, when the hysteretic damping is null, a viscous damping component is always considered.This overcomes the limitations of other computer codes that result in unrealistic amplifications in the case of low energy earthquakes [32,41].The Rayleigh damping target frequencies associated to the target soil damping at low strains, used in this work, are those commonly recommended in DEEPSOIL: (a) the natural frequency of the soil profile (0.69 Hz) and (b) 5 times the natural frequency of the soil profile (3.46 Hz).
In Figures 10 and 11 the comparison between the experimental and the adopted normalized modulus reduction and damping ratio curves corresponding to the A1 and Upper Pancone (B1-B5) layers are presented.The damping ratio curves shown in Figures 10 and 11 related with ONDA don't include the viscous damping at low strain levels that has been obviously included in the analyses to avoid zero damping values.ONDA overestimates the hysteretic damping at high strain levels compared to DEEPSOIL.This fact justifies the slightly lower spectral acceleration values obtained with ONDA and shown in the following section.
a linear combination of the mass and the stiffness matrixes (i.e., Rayleigh damping formulation).Consequently, even at small strain levels, when the hysteretic damping is null, a viscous damping component is always considered.This overcomes the limitations of other computer codes that result in unrealistic amplifications in the case of low energy earthquakes [32,41].The Rayleigh damping target frequencies associated to the target soil damping at low strains, used in this work, are those commonly recommended in DEEPSOIL: (a) the natural frequency of the soil profile (0.69 Hz) and (b) 5 times the natural frequency of the soil profile (3.46 Hz).
In Figures 10 and 11 the comparison between the experimental and the adopted normalized modulus reduction and damping ratio curves corresponding to the A1 and Upper Pancone (B1-B5) layers are presented.The damping ratio curves shown in Figures 10 and 11 related with ONDA don't include the viscous damping at low strain levels that has been obviously included in the analyses to avoid zero damping values.ONDA overestimates the hysteretic damping at high strain levels compared to DEEPSOIL.This fact justifies the slightly lower spectral acceleration values obtained with ONDA and shown in the following section.a linear combination of the mass and the stiffness matrixes (i.e., Rayleigh damping formulation).
Consequently, even at small strain levels, when the hysteretic damping is null, a viscous damping component is always considered.This overcomes the limitations of other computer codes that result in unrealistic amplifications in the case of low energy earthquakes [32,41].The Rayleigh damping target frequencies associated to the target soil damping at low strains, used in this work, are those commonly recommended in DEEPSOIL: (a) the natural frequency of the soil profile (0.69 Hz) and (b) 5 times the natural frequency of the soil profile (3.46 Hz).
In Figures 10 and 11 the comparison between the experimental and the adopted normalized modulus reduction and damping ratio curves corresponding to the A1 and Upper Pancone (B1-B5) layers are presented.The damping ratio curves shown in Figures 10 and 11 related with ONDA don't include the viscous damping at low strain levels that has been obviously included in the analyses to avoid zero damping values.ONDA overestimates the hysteretic damping at high strain levels compared to DEEPSOIL.This fact justifies the slightly lower spectral acceleration values obtained with ONDA and shown in the following section.

Seismic Response Analysis Results
Figure 12 shows the difference between the soil shear strength profile obtained with the commonly used MKZ model and that computed with Equation ( 8), used to perform the nonlinear analyses with strength-controlled constitutive models (DEEPSOIL GQ/H and ONDA).

Seismic Response Analysis results
Figure 12 shows the difference between the soil shear strength profile obtained with the commonly used MKZ model and that computed with Equation ( 8), used to perform the nonlinear analyses with strength-controlled constitutive models (DEEPSOIL GQ/H and ONDA).The MKZ model is not able to correctly reproduce the shear strength profile and largely overestimate the soil resistance at failure down to 4 m in depth, while underestimate the soil resistance at higher depths.This affects the nonlinear soil response analysis results.
The results for a return period of 130 and 500 years are reported in Figures 13 and 14 in terms of horizontal acceleration response spectra computed at ground surface (averaged over all the input motions applied).Average input spectrum at seismic bedrock and amplified spectra are also reported.The MKZ model is not able to correctly reproduce the shear strength profile and largely overestimate the soil resistance at failure down to 4 m in depth, while underestimate the soil resistance at higher depths.This affects the nonlinear soil response analysis results.
The results for a return period of 130 and 500 years are reported in Figures 13 and 14 in terms of horizontal acceleration response spectra computed at ground surface (averaged over all the input motions applied).Average input spectrum at seismic bedrock and amplified spectra are also reported.

Seismic Response Analysis results
Figure 12 shows the difference between the soil shear strength profile obtained with the commonly used MKZ model and that computed with Equation ( 8), used to perform the nonlinear analyses with strength-controlled constitutive models (DEEPSOIL GQ/H and ONDA).The MKZ model is not able to correctly reproduce the shear strength profile and largely overestimate the soil resistance at failure down to 4 m in depth, while underestimate the soil resistance at higher depths.This affects the nonlinear soil response analysis results.
The results for a return period of 130 and 500 years are reported in Figures 13 and 14 in terms of horizontal acceleration response spectra computed at ground surface (averaged over all the input motions applied).Average input spectrum at seismic bedrock and amplified spectra are also reported.
(a) (b) Figure 15a,b show the Fourier spectrum ratios between the ground surface and the base of the soil profile for some representative accelerograms at return periods of 130 and 500 years, respectively.In both the Figures and for all the accelerograms, the ratio is generally less than 1 for frequencies above 2 Hz in case of nonlinear analysis.Figure 15a,b show the Fourier spectrum ratios between the ground surface and the base of the soil profile for some representative accelerograms at return periods of 130 and 500 years, respectively.In both the Figures and for all the accelerograms, the ratio is generally less than 1 for frequencies above 2 Hz in case of nonlinear analysis.Figure 16 shows the comparisons between the average response spectra, while Figure 17 the average surface-to-bedrock amplification function for both the return periods considered in this study.The SRA comparisons shown in Figures 16 and 17 reveal the importance of the selected analysis method (equivalent-linear or nonlinear) in the definition of pseudo-acceleration.
The comparison between equivalent-linear and nonlinear results revealed more evident differences.In fact, equivalent-linear method overestimates the ordinates between 0 and 1.0 s period for 130 years RP, whereas for 500 years RP ordinates are overestimated between 0 and 1.5 s.Maximum overestimation is above 80% in the interval 0-0.5 s.
The equivalent-linear models (e.g., STRATA [28]) consider a single elastic modulus referred to a representative strain level (65% or M-1/10 of the max range).
In nonlinear analyses, on the other hand, the elastic modulus G is updated at each time-step and is reduced considerably for medium-large deformations, causing an increase in the deposit fundamental period and thus shifting the peaks of the spectrum to higher periods.
In the case of 130 years RP above 0.6 s, the nonlinear analyses practically coincide despite the fact that the formulation of the backbone curve adopted in ONDA [32,33] (Ramberg-Osgood) and DEEPSOIL [29] (MKZ [8] and GQ/H [30,31]) is different.
In the case of 500 years RP the differences between the commonly used nonlinear MKZ model and strength-controlled constitutive models (GQ/H and ONDA) become higher compared to seismic motions for 130 years RP.The MKZ model overestimates the ordinates obtained with ONDA of about 10-20% between 0 and 0.8 s.This is due to the different soil modelling for medium-large shear strain levels.
The latter aspect can be clarified showing in Figure 18 the peak shear strain profiles for 130 and 500 years RP seismic motions.In case of 500 years RP acceleration time histories, the attained shear The SRA comparisons shown in Figures 16 and 17 reveal the importance of the selected analysis method (equivalent-linear or nonlinear) in the definition of pseudo-acceleration.
The comparison between equivalent-linear and nonlinear results revealed more evident differences.In fact, equivalent-linear method overestimates the ordinates between 0 and 1.0 s period for 130 years RP, whereas for 500 years RP ordinates are overestimated between 0 and 1.5 s.Maximum overestimation is above 80% in the interval 0-0.5 s.
The equivalent-linear models (e.g., STRATA [28]) consider a single elastic modulus referred to a representative strain level (65% or M-1/10 of the max range).
In nonlinear analyses, on the other hand, the elastic modulus G is updated at each time-step and is reduced considerably for medium-large deformations, causing an increase in the deposit fundamental period and thus shifting the peaks of the spectrum to higher periods.
In the case of 130 years RP above 0.6 s, the nonlinear analyses practically coincide despite the fact that the formulation of the backbone curve adopted in ONDA [32,33] (Ramberg-Osgood) and DEEPSOIL [29] (MKZ [8] and GQ/H [30,31]) is different.
In the case of 500 years RP the differences between the commonly used nonlinear MKZ model and strength-controlled constitutive models (GQ/H and ONDA) become higher compared to seismic motions for 130 years RP.The MKZ model overestimates the ordinates obtained with ONDA of about 10-20% between 0 and 0.8 s.This is due to the different soil modelling for medium-large shear strain levels.
The latter aspect can be clarified showing in Figure 18 the peak shear strain profiles for 130 and 500 years RP seismic motions.In case of 500 years RP acceleration time histories, the attained shear The SRA comparisons shown in Figures 16 and 17 reveal the importance of the selected analysis method (equivalent-linear or nonlinear) in the definition of pseudo-acceleration.
The comparison between equivalent-linear and nonlinear results revealed more evident differences.In fact, equivalent-linear method overestimates the ordinates between 0 and 1.0 s period for 130 years RP, whereas for 500 years RP ordinates are overestimated between 0 and 1.5 s.Maximum overestimation is above 80% in the interval 0-0.5 s.
The equivalent-linear models (e.g., STRATA [28]) consider a single elastic modulus referred to a representative strain level (65% or M-1/10 of the max range).
In nonlinear analyses, on the other hand, the elastic modulus G is updated at each time-step and is reduced considerably for medium-large deformations, causing an increase in the deposit fundamental period and thus shifting the peaks of the spectrum to higher periods.
In the case of 130 years RP above 0.6 s, the nonlinear analyses practically coincide despite the fact that the formulation of the backbone curve adopted in ONDA [32,33] (Ramberg-Osgood) and DEEPSOIL [29] (MKZ [8] and GQ/H [30,31]) is different.
In the case of 500 years RP the differences between the commonly used nonlinear MKZ model and strength-controlled constitutive models (GQ/H and ONDA) become higher compared to seismic motions for 130 years RP.The MKZ model overestimates the ordinates obtained with ONDA of about 10-20% between 0 and 0.8 s.This is due to the different soil modelling for medium-large shear strain levels.
The latter aspect can be clarified showing in Figure 18 the peak shear strain profiles for 130 and 500 years RP seismic motions.In case of 500 years RP acceleration time histories, the attained shear strain levels are larger compared to 130 years RP input motions.For such a reason becomes relevant to consider a strength-controlled constitutive model for the soil response.

Closing Remarks
The paper shows the results of free-field seismic response analyses (SRA) carried out considering the subsoil of Piazza dei Miracoli in Pisa, where the Leaning Tower is located.Given that this site has been subject to a considerable number of investigations, there is a large number of data available for soil characterization.
The shear wave velocity profile has been recently extended down to more than 100 m below the ground level [6].One-dimensional SRA were carried out by means of three computer codes: STRATA [28], DEEPSOIL [29] and ONDA [32,33].
The first code performs the analysis in the frequency domain considering an equivalent-linear soil model, whereas DEEPSOIL and ONDA analyze the problem in the time domain assuming true nonlinear soil behavior.All the parameters of the models were calibrated using the same soil data, therefore all three models can be considered as being equivalent.
As discussed in [6], seismic input is affected by some dispersion due to uncertainties related to lack of knowledge.The SRA analyses are also affected by aleatory uncertainty.
The main target of this paper is to show the influence of using new strength-controlled constitutive models in case of 1D nonlinear response analyses.In fact, most of computer codes currently available for nonlinear seismic response analysis use constitutive models able to capture small-strain behavior, but the large-strain shear strength is left uncontrolled.This can affect the assessment of a 1-D response analysis, as in the Pisa subsoil conditions.Nonlinear models seem more appropriate in reproducing the true behavior of the subsoil, however calibrating the parameters seems more difficult.This highlights the importance of considering a more proper and realistic soil strength profile in nonlinear analysis methods starting from medium-severe seismic actions.

Closing Remarks
The paper shows the results of free-field seismic response analyses (SRA) carried out considering the subsoil of Piazza dei Miracoli in Pisa, where the Leaning Tower is located.Given that this site has been subject to a considerable number of investigations, there is a large number of data available for soil characterization.
The shear wave velocity profile has been recently extended down to more than 100 m below the ground level [6].One-dimensional SRA were carried out by means of three computer codes: STRATA [28], DEEPSOIL [29] and ONDA [32,33].
The first code performs the analysis in the frequency domain considering an equivalent-linear soil model, whereas DEEPSOIL and ONDA analyze the problem in the time domain assuming true nonlinear soil behavior.All the parameters of the models were calibrated using the same soil data, therefore all three models can be considered as being equivalent.
As discussed in [6], seismic input is affected by some dispersion due to uncertainties related to lack of knowledge.The SRA analyses are also affected by aleatory uncertainty.
The main target of this paper is to show the influence of using new strength-controlled constitutive models in case of 1D nonlinear response analyses.In fact, most of computer codes currently available for nonlinear seismic response analysis use constitutive models able to capture small-strain behavior, but the large-strain shear strength is left uncontrolled.This can affect the assessment of a 1-D response analysis, as in the Pisa subsoil conditions.Nonlinear models seem more appropriate in reproducing the true behavior of the subsoil, however calibrating the parameters seems more difficult.
In this work, differences between the commonly used nonlinear MKZ model and strength-controlled constitutive models (GQ/H and the model implemented in ONDA) are presented.The MKZ model overestimates the pseudo-acceleration values obtained with ONDA and the GQ/H model, especially in case of 500 years RP seismic input motions.This is due to the different soil modelling for medium-large shear strain levels.In fact, the attained shear strain levels are larger compared to 130 years RP acceleration time histories.This reveals the importance to take into account of a proper soil strength profile in nonlinear analysis methods starting from medium-severe seismic actions.Higher differences were found between equivalent-linear and nonlinear codes.The maximum differences between equivalent-linear and nonlinear models are above 80% in both the 130 and 500 years return periods.For RP = 130 years, the differences decrease and become negligible starting from a period approximately equal to 1 s, whereas the differences remain around 40% for RP = 500 years.
The differences among the mean spectra obtained by means of different soil response modelling (equivalent-linear, MKZ model, GQ/H model and ONDA) might be not relevant for a structural designer in case of a new standard building, since it is possible to employ code-defined spectrum or, as an alternative, select the most conservative one.
The selection of the appropriate spectrum assumes particular relevance in case of historical buildings or monuments, because in this case structural interventions have to be respectful of construction integrity and a selection of excessively conservative actions could lead to inappropriate operations [17,42].Moreover, in the case of the Leaning Tower of Pisa, since the fundamental period is about one second, the uncertainties related to nonlinear models seem negligible with respect to the uncertainty related to seismic input.Further studies and investigations have been planned in order to improve the knowledge of subsoil below 100 m depth.This is expected to improve the capabilities of considered models to forecast the ground response.

Figure 2 .
Figure 2. Morphological map of the study area.The scale is the same as the previous figure.(Adapted from: Tuscany Region Database-Geoscopio [10]).

Figure 2 .
Figure 2. Morphological map of the study area.The scale is the same as the previous figure.(Adapted from: Tuscany Region Database-Geoscopio [10]).

Figure 2 .
Figure 2. Morphological map of the study area.The scale is the same as the previous figure.(Adapted from: Tuscany Region Database-Geoscopio [10]).
Horizon B consists primarily of marine clay, which extends down to a depth of about 40 m.It is subdivided into four distinct layers.The upper layer is soft sensitive clay, locally known as Pancone.It is underlain by an intermediate layer of stiffer clay, which in turn overlies a sand layer (the intermediate sand).The bottom layer of horizon B is a normally consolidated clay known as the lower clay.Horizon B is very uniform laterally near the tower.Geosciences 2018, 8, x FOR PEER REVIEW 4 of 20
and an SDMT test (Seismic Dilatometer Test).The results obtained from these tests show a satisfactory agreement of the Vs values.All the tests reached a maximum depth of 40 m except for the CH test, which reached a depth of 65 m.None of the tests was successful in identifying the seismic bedrock.
and an SDMT test (Seismic Dilatometer Test).The results obtained from these tests show a satisfactory agreement of the Vs values.All the tests reached a maximum depth of 40 m except for the CH test, which reached a depth of 65 m.None of the tests was successful in identifying the seismic bedrock.

Figure 5 .
Figure 5. Profile of the shear wave velocity.

Figure 5 .
Figure 5. Profile of the shear wave velocity.

Figure 6 .
Figure 6.Time histories of the scaled accelerograms for an RP of 130 years.Figure 6.Time histories of the scaled accelerograms for an RP of 130 years.

Figure 6 .
Figure 6.Time histories of the scaled accelerograms for an RP of 130 years.Figure 6.Time histories of the scaled accelerograms for an RP of 130 years.

Figure 7 .
Figure 7. Time histories of the scaled accelerograms for an RP of 500 years.

Figure 8 Figure 8 .
Figure8reports the selected input motion spectra and the results in terms of average response spectrum for 130 years (left) and 500 years (right).The dispersion of the records (average spectrum (μ) ± 1 standard deviation (σ)) is reported for the two cases.For RP = 130 years, the dispersion is large, with a maximum COV (σ/μ) of 0.996 in the range 0.31-1 s, while for RP = 500 years there is a large dispersion too, but slightly lower compared to RP = 130 years, and a maximum COV of 0.643 in the above cited range.

Figure 7 .
Figure 7. Time histories of the scaled accelerograms for an RP of 500 years.

Figure 8 20 Figure 7 .
Figure8reports the selected input motion spectra and the results in terms of average response spectrum for 130 years (left) and 500 years (right).The dispersion of the records (average spectrum (µ) ± 1 standard deviation (σ)) is reported for the two cases.For RP = 130 years, the dispersion is large, with a maximum COV (σ/µ) of 0.996 in the range 0.31-1 s, while for RP = 500 years there is a large dispersion too, but slightly lower compared to RP = 130 years, and a maximum COV of 0.643 in the above cited range.

Figure 8 Figure 8 .
Figure8reports the selected input motion spectra and the results in terms of average response spectrum for 130 years (left) and 500 years (right).The dispersion of the records (average spectrum (μ) ± 1 standard deviation (σ)) is reported for the two cases.For RP = 130 years, the dispersion is large, with a maximum COV (σ/μ) of 0.996 in the range 0.31-1 s, while for RP = 500 years there is a large dispersion too, but slightly lower compared to RP = 130 years, and a maximum COV of 0.643 in the above cited range.

Figure 11 .
Figure 11.Comparison between the experimental and the adopted: (a) modulus reduction curve; (b) damping ratio curve (Upper Pancone).In this Figure, triangles represent the best-fit of the collected experimental data related to the Upper Pancone as reported in [7].

Figure 11 .
Figure 11.Comparison between the experimental and the adopted: (a) modulus reduction curve; (b) damping ratio curve (Upper Pancone).In this Figure, triangles represent the best-fit of the collected experimental data related to the Upper Pancone as reported in [7].

Figure 11 .
Figure 11.Comparison between the experimental and the adopted: (a) modulus reduction curve; (b) damping ratio curve (Upper Pancone).In this Figure, triangles represent the best-fit of the collected experimental data related to the Upper Pancone as reported in [7].

Figure
Figure15a,b show the Fourier spectrum ratios between the ground surface and the base of the soil profile for some representative accelerograms at return periods of 130 and 500 years, respectively.In both the Figures and for all the accelerograms, the ratio is generally less than 1 for frequencies above 2 Hz in case of nonlinear analysis.

Figure 15 .
Figure 15.Fourier spectrum ratios between the surface and the base of the soil profile for some relevant input motions: (a) return period = 130 years; (b) return period = 500 years.

Figure 15 .
Figure 15.Fourier spectrum ratios between the surface and the base of the soil profile for some relevant input motions: (a) return period = 130 years; (b) return period = 500 years.

Figure 16
Figure16shows the comparisons between the average response spectra, while Figure17the average surface-to-bedrock amplification function for both the return periods considered in this study.

Table 1 .
Subsoil model adopted for site response analysis. Hor.

Table 1 .
Subsoil model adopted for site response analysis.