Effects of Soil-Foundation-Interaction on the Seismic Response of a Cooling Tower by 3D-FEM Analysis

This paper reports on the results of soil-foundation numerical modelling and the seismic response of a cooling tower founded on piles of a petrochemical facility located in the city of Augusta (Sicily, Italy). The city was affected in the past by some destructive earthquakes (1693, 1848, and 1990) that damaged a large territory of Southeastern Sicily, which was characterized by a very high seismic hazard. With this aim, the paper reports the FEM modelling of the seismic behaviour of the coupled soil-structure system. To determine the soil profile and the geotechnical characteristics, laboratory and in situ investigations were carried out in the studied area. The seismic event occurred in January 1693 and has been chosen as a scenario earthquake. Moreover, a parametric study with different input motions has also been carried out. A Mohr-Coulomb model has been adopted for the soil, and structural elements have been simulated by means of an elastic constitutive model. Two different vertical alignments have been analysed, considering both the free-field condition and the soil-structure interaction. The dynamic response has been investigated in terms of accelerations, response spectra, and amplification functions. The results have also been compared with those provided by Italian technical regulations. Finally, the seismic response of the coupled soil-structure system has been further examined in terms of peak bending moments along the pile foundation, emphasizing the possibility of a kinematic interaction on piles induced by the seismic action.

In an industrial plant, a seismic event can cause not only economic impacts due to the interruption of the production cycle, but also casualties as a consequence of components collapsing, explosions, fires, and a release of toxic substances [1,2]. In addition, the collapse of an apparatus can induce the simultaneous damage of other apparatuses, triggering an uncontrolled domino effect [1,3].
To evaluate the seismic vulnerability of petrochemical facilities, several codes and guidelines are used. Many of these codes have been developed primarily for building design and generally offer insufficiently detailed guidelines for the design and the evaluation of structures commonly found in petrochemical facilities. Consequently, many petrochemical operating companies have developed their own internal standards and guidelines for addressing the seismic design and evaluation issues.
In Italy, about one-third of relevant risk plants are located in medium to high seismic areas [3]. For this reason, several research studies were carried out with the aim to analyze the seismic responses of components in industrial process plants [4]. The city has been destroyed by several disastrous earthquakes with an MCS intensity from VI to X ( Figure 2).
Geosciences 2020, 10, x FOR PEER REVIEW 4 50 and 300 m. The upper part consists of recent alluvial soil that overhangs Pleisto marine clay, called Augusta blue-grey clay [15]. Figure 1 shows the geological sketch of the Hyblean Plateau. The city has been destroyed by several disastrous earthquakes with an MCS inten from VI to X ( Figure 2). The earthquake on 11 January 1693, known as "Val di Noto earthquake," stru vast territory of Southeastern Sicily. It had an estimated magnitude of 7.4 on the mom magnitude scale and caused the partial, and, in many cases, total, destruction of 57 c and 60,000 casualties [18,19,20].
The Etna earthquake that took place on 20 February 1818 caused moderate dam The earthquake on 11 January 1693, known as "Val di Noto earthquake," struck a vast territory of Southeastern Sicily. It had an estimated magnitude of 7.4 on the moment magnitude scale and caused the partial, and, in many cases, total, destruction of 57 cities and 60,000 casualties [18][19][20].
The Etna earthquake that took place on 20 February 1818 caused moderate damage, but its effects were noticed over a vast area [21]. Some localities of Southeastern Sicily were shaken by the 1848 earthquake. The largest damage was observed in Augusta, where almost two-thirds of the buildings collapsed [22].
On December 1908, a devastating earthquake occurred along the Strait of Messina between the eastern tip of Sicily and the western tip of Calabria in Southern Italy. The greatest damage was experienced near the Strait of Messina. In Augusta, the violence reached VI-VII degrees on the Mercalli scale.
In the studied area, a large number of buildings exhibited fierce damage or collapsed due to the earthquake of 13 December 1990. Although the shock was of a moderate magnitude (M W = 5.4), it caused deaths, injuries, and widespread damage due to the soft-soil conditions and construction deficiencies [23].
The seismic event that occurred in January 1693 has been chosen as the scenario earthquake. Figure 3 shows the intensity map for the 1693 earthquake.
Geosciences 2020, 10, x FOR PEER REVIEW 5 of 32 In the studied area, a large number of buildings exhibited fierce damage or collapsed due to the earthquake of 13 December 1990. Although the shock was of a moderate magnitude (MW = 5.4), it caused deaths, injuries, and widespread damage due to the softsoil conditions and construction deficiencies [23].
The seismic event that occurred in January 1693 has been chosen as the scenario earthquake. Figure 3 shows the intensity map for the 1693 earthquake.  [24,25], modified, INGV database).

Geotechnical Soil Properties
The industrial district of Augusta is located on a flat area of the coastal hinterland. The area faces the Gulf of Augusta and consists of production facilities, including mainly petrochemical and cement factories ( Figure 4).  [24,25], modified, INGV database).

Geotechnical Soil Properties
The industrial district of Augusta is located on a flat area of the coastal hinterland. The area faces the Gulf of Augusta and consists of production facilities, including mainly petrochemical and cement factories ( Figure 4).  [24,25], modified, INGV database).

Geotechnical Soil Properties
The industrial district of Augusta is located on a flat area of the coastal hinterland. The area faces the Gulf of Augusta and consists of production facilities, including mainly petrochemical and cement factories ( Figure 4).    The locations of the boreholes, driven to a depth of 30 m, are reported in Figure 5.  Figure 5. The results of the boreholes are as follows.  The results of the boreholes are as follows. The water table lies around 11 m below the ground surface. The boreholes were equipped for the execution of the Down Hole tests. In Figure 6, the shear and compression wave velocities are shown against the depth, for borehole S4, as an example.  The investigation reached a depth of 30 m. The VS values ranged from 160 to 360 m/s. Higher values were reported for depths between 11 and 19 m.
The dynamic Poisson ratio variation (ν) with depth, obtained from the DH test (borehole S4), is reported in Figure 7. The ν values range from 0.417 to 0.496. The investigation reached a depth of 30 m. The V S values ranged from 160 to 360 m/s. Higher values were reported for depths between 11 and 19 m.
The dynamic Poisson ratio variation (ν) with depth, obtained from the DH test (borehole S4), is reported in Figure 7. The ν values range from 0.417 to 0.496. In addition to in situ investigations, laboratory tests were carried out on undisturbed samples, as reported in Table 1. Based on the laboratory tests, index properties and strength parameters of the deposits are reported in Table 2.  In addition to in situ investigations, laboratory tests were carried out on undisturbed samples, as reported in Table 1. Based on the laboratory tests, index properties and strength parameters of the deposits are reported in Table 2.
Shear modulus G and damping ratio D were obtained from resonant column tests (RCT) and cyclic loading torsional shear tests (CLTST) performed by means of a resonant column/cyclic loading torsional shear apparatus. The same specimen was first subject to RCT, and then to CLTST, after a rest period of 24 h with opened drainage [26][27][28][29].
Values  Shear modulus G and damping ratio D were obtained from resonant column tests (RCT) and cyclic loading torsional shear tests (CLTST) performed by means of a resonant column/cyclic loading torsional shear apparatus. The same specimen was first subject to RCT, and then to CLTST, after a rest period of 24 h with opened drainage [26,27,28,29].
Values   According to these data, it is possible to observe that higher values of D are obtained by CLTS than RCT.

Petrochemical Facilities
Southeastern Sicily is considered one of the zones of Italy with the highest seismic risk. In the city of Augusta, historically affected by strong earthquakes and high-intensity tsunami events, there are important industrial petrochemical facilities. They are classified with high risk, according to the Seveso III Directive on the control of major-accident hazards involving dangerous substances [30].
One of the main changes introduced by the Seveso III Directive is the recognition of the impact of natural risks on safety management in major accident-generated industrial plants. The accidents generated are called NaTech (Natural Hazard Triggering According to these data, it is possible to observe that higher values of D are obtained by CLTS than RCT.

Petrochemical Facilities
Southeastern Sicily is considered one of the zones of Italy with the highest seismic risk. In the city of Augusta, historically affected by strong earthquakes and high-intensity tsunami events, there are important industrial petrochemical facilities. They are classified with high risk, according to the Seveso III Directive on the control of major-accident hazards involving dangerous substances [30].
One of the main changes introduced by the Seveso III Directive is the recognition of the impact of natural risks on safety management in major accident-generated industrial plants. The accidents generated are called NaTech (Natural Hazard Triggering Technological Disasters). The natural events can be rainfall, floods, earthquake, volcanic phenomena, or hurricane [31]. Earthquakes can cause extensive damage to industrial facilities [32].
In this paper, the effects of the soil-foundation-interaction on the seismic response of a cooling tower, located in Augusta, have been studied by 3D-FEM analysis.
The cooling tower under consideration, called C-352, consists of the following components: reinforced concrete (RC) piles, RC plate, RC pedestal, and equipment. Figure 11 shows a scheme of the cooling tower. phenomena, or hurricane [31]. Earthquakes can cause extensive damage to industrial facilities [32]. In this paper, the effects of the soil-foundation-interaction on the seismic response of a cooling tower, located in Augusta, have been studied by 3D-FEM analysis.
The cooling tower under consideration, called C-352, consists of the following components: reinforced concrete (RC) piles, RC plate, RC pedestal, and equipment. Figure  11 shows a scheme of the cooling tower.

Model Geometry and Boundary Conditions of Soil
From in situ characterization (boreholes S1, S2, S3, and S4), it has been possible to define the following soil stratigraphy: Backfill (from 0.00 to 2.50 m), Alluvial Deposits (from 2.50 to 12.00 m), Sandy Clay (from 12.00 to 17.00 m), and Blue-Grey Clay (from 17.00 to 30.00 m). The soil profile is shown in Figure 14.

Model Geometry and Boundary Conditions of Soil
From in situ characterization (boreholes S1, S2, S3, and S4), it has been possible to define the following soil stratigraphy: Backfill (from 0.00 to 2.50 m), Alluvial Deposits (from 2.50 to 12.00 m), Sandy Clay (from 12.00 to 17.00 m), and Blue-Grey Clay (from 17.00 to 30.00 m). The soil profile is shown in Figure 14.

Model Geometry and Boundary Conditions of Soil
From in situ characterization (boreholes S1, S2, S3, and S4), it has been possible to define the following soil stratigraphy: Backfill (from 0.00 to 2.50 m), Alluvial Deposits (from 2.50 to 12.00 m), Sandy Clay (from 12.00 to 17.00 m), and Blue-Grey Clay (from 17.00 to 30.00 m). The soil profile is shown in Figure 14. The dynamic response model requires knowledge of bedrock depth. Marzorati et al. [33] and Fiorini et al. [34] carried out a noise measurements campaign in order to obtain the soil resonant frequency in the industrial area. The measures provided the resonant frequency of the soil deposit (f0) and gave indications about the average shear wave velocity and the depth of seismic bedrock. The analysis of the noise measurements using the spectral ratios HVSR (Horizontal to Vertical Spectral Ratio) showed the lowest values of f0 (0.69-1.37 Hz) along the coastal area under consideration ( Figure 5). Moreover, the results indicated that the most important geotechnical unit in terms of soil response is the blue clay with an average shear wave velocity of 600 m/s. The same value of VS is also reported by Tortorici [35]. Finally, the results allowed us to identify the surface that The dynamic response model requires knowledge of bedrock depth. Marzorati et al. [33] and Fiorini et al. [34] carried out a noise measurements campaign in order to obtain the soil resonant frequency in the industrial area. The measures provided the resonant frequency of the soil deposit (f 0 ) and gave indications about the average shear wave velocity and the depth of seismic bedrock. The analysis of the noise measurements using the spectral ratios HVSR (Horizontal to Vertical Spectral Ratio) showed the lowest values of f 0 (0.69-1.37 Hz) along the coastal area under consideration ( Figure 5). Moreover, the results indicated that the most important geotechnical unit in terms of soil response is the blue clay with an average shear wave velocity of 600 m/s. The same value of V S is also reported by Tortorici [35]. Finally, the results allowed us to identify the surface that generates the amplifications at the base of the blue-grey clays, which constitutes the seismic bedrock [33,34]. The thickness of the blue-grey clay is equal to 73 m in the studied area [36]. The soil model is shown in Table 3, where V S values are obtained from D-H tests and from previous considerations. The values of the unit weight have been provided by the laboratory tests. Finally, the water head has been imposed at a depth of 11 m. The fundamental natural frequency of the soil deposit can be calculated according to the following formula.
where V S is the shear wave velocity in the soil deposit of thickness H. The fundamental natural frequency of the soil deposit is equal to 1.37 Hz.
The dynamic input has been applied along the x-direction by imposing a given displacement at the bottom boundary. In order to minimize boundary effects as much as possible, the length of the model L (along the x-direction) has been chosen to equal 180 m. Instead, the boundary efreects are minimal along the y-direction. Therefore, the width of the model W has been chosen to equal nine times the equipment's footprint, i.e., 31.5 m.
The basic soil elements of the 3D finite element mesh are the 10-node tetrahedral elements. The average element size (AES) respects the following equation [37]. AES ≤ V s,min 6 ÷ 8f max (2) where V S,min is the lowest wave velocity and f max is the maximum frequency component of the input wave. The number of elements is 308,703. The finite element mesh is shown in Figure 15.
Geosciences 2020, 10, x FOR PEER REVIEW 14 of 32 generates the amplifications at the base of the blue-grey clays, which constitutes the seismic bedrock [33,34]. The thickness of the blue-grey clay is equal to 73 m in the studied area [36]. The soil model is shown in Table 3, where VS values are obtained from D-H tests and from previous considerations. The values of the unit weight have been provided by the laboratory tests. Finally, the water head has been imposed at a depth of 11 m. The fundamental natural frequency of the soil deposit can be calculated according to the following formula.
where VS is the shear wave velocity in the soil deposit of thickness H. The fundamental natural frequency of the soil deposit is equal to 1.37 Hz. The dynamic input has been applied along the x-direction by imposing a given displacement at the bottom boundary. In order to minimize boundary effects as much as possible, the length of the model L (along the x-direction) has been chosen to equal 180 m. Instead, the boundary efreects are minimal along the y-direction. Therefore, the width of the model W has been chosen to equal nine times the equipment's footprint, i.e., 31.5 m.
The basic soil elements of the 3D finite element mesh are the 10-node tetrahedral elements. The average element size (AES) respects the following equation [37].
AES ≤ V , 6 ÷ 8f (2) where VS,min is the lowest wave velocity and fmax is the maximum frequency component of the input wave. The number of elements is 308,703.
The finite element mesh is shown in Figure 15. With regard to the boundary conditions, to generate the initial stresses due to gravity loading, default fixities have been applied, corresponding to normally fixed vertical boundaries and a fully fixed base. For the dynamic phase, the free-field and compliant base boundary conditions have been considered ( Figure 16). With regard to the boundary conditions, to generate the initial stresses due to gravity loading, default fixities have been applied, corresponding to normally fixed vertical boundaries and a fully fixed base. For the dynamic phase, the free-field and compliant base boundary conditions have been considered ( Figure 16). The free field condition simulates the propagation of waves into the far field with a minimum reflection at the boundary. The domain is reduced to the area of interest and the free field motion is applied to the boundaries employing free-field elements. The motion is transferred from the free field elements to the main domain by applying equivalent normal and shear forces. Two dashpots are added in the normal and shear direction at each node of the model boundary ( Figure 17) [38,39]. The compliant base boundary [40] simulates the continuation of waves into the deep soil with a minimum reflection at the bottom boundary. Free field elements can be attached to the bottom of the main domain and the input signal is transferred to the main domain by applying equivalent normal and shear forces [38,39].
Model geometry and boundary conditions are shown in Figure 18. The free field condition simulates the propagation of waves into the far field with a minimum reflection at the boundary. The domain is reduced to the area of interest and the free field motion is applied to the boundaries employing free-field elements. The motion is transferred from the free field elements to the main domain by applying equivalent normal and shear forces. Two dashpots are added in the normal and shear direction at each node of the model boundary ( Figure 17) [38,39]. The free field condition simulates the propagation of waves into the far field with a minimum reflection at the boundary. The domain is reduced to the area of interest and the free field motion is applied to the boundaries employing free-field elements. The motion is transferred from the free field elements to the main domain by applying equivalent normal and shear forces. Two dashpots are added in the normal and shear direction at each node of the model boundary ( Figure 17) [38,39]. The compliant base boundary [40] simulates the continuation of waves into the deep soil with a minimum reflection at the bottom boundary. Free field elements can be attached to the bottom of the main domain and the input signal is transferred to the main domain by applying equivalent normal and shear forces [38,39].
Model geometry and boundary conditions are shown in Figure 18. The compliant base boundary [40] simulates the continuation of waves into the deep soil with a minimum reflection at the bottom boundary. Free field elements can be attached to the bottom of the main domain and the input signal is transferred to the main domain by applying equivalent normal and shear forces [38,39].
Model geometry and boundary conditions are shown in Figure 18.
Geosciences 2020, 10, x FOR PEER REVIEW 16 of Figure 18. Model geometry and boundary conditions.

Full-Coupled FEM Model
Dynamic response of a structure deviates from the free field (FF) site respon analysis because of a kinematic and inertial interaction [11,41]. The dynamic response a fully-coupled soil-structure system has been performed by means of PLAXIS software. Two different vertical alignments have been analyzed: along the structure (S and far from the structure (FF), as shown in Figure 19.

Mohr-Coulomb Model
The Mohr-Coulomb model has been used for both static and dynamic calculatio The linear elastic perfectly-plastic model requires five parameters as input: Youn modulus (E) and Poisson's ratio (ν) for soil elasticity, friction angle (φ') and cohesion for soil plasticity, and the dilatancy angle (ψ). The stiffness has been defined based wave velocity obtained from the geophysical tests. The following relationships have be used to determine the Young's modulus (E) and Poisson's ratio (ν) from the wa velocities.

Full-Coupled FEM Model
Dynamic response of a structure deviates from the free field (FF) site response analysis because of a kinematic and inertial interaction [11,41]. The dynamic response of a fullycoupled soil-structure system has been performed by means of PLAXIS 3D software. Two different vertical alignments have been analyzed: along the structure (SSI) and far from the structure (FF), as shown in Figure 19.

Full-Coupled FEM Model
Dynamic response of a structure deviates from the free field (FF) site response analysis because of a kinematic and inertial interaction [11,41]. The dynamic response of a fully-coupled soil-structure system has been performed by means of PLAXIS 3D software. Two different vertical alignments have been analyzed: along the structure (SSI) and far from the structure (FF), as shown in Figure 19.

Mohr-Coulomb Model
The Mohr-Coulomb model has been used for both static and dynamic calculations. The linear elastic perfectly-plastic model requires five parameters as input: Young's modulus (E) and Poisson's ratio (ν) for soil elasticity, friction angle (φ') and cohesion (c') for soil plasticity, and the dilatancy angle (ψ). The stiffness has been defined based on wave velocity obtained from the geophysical tests. The following relationships have been used to determine the Young's modulus (E) and Poisson's ratio (ν) from the wave velocities.

Mohr-Coulomb Model
The Mohr-Coulomb model has been used for both static and dynamic calculations. The linear elastic perfectly-plastic model requires five parameters as input: Young's modulus (E) and Poisson's ratio (ν) for soil elasticity, friction angle (ϕ') and cohesion (c') for soil plasticity, and the dilatancy angle (ψ). The stiffness has been defined based on wave velocity obtained from the geophysical tests. The following relationships have been used to determine the Young's modulus (E) and Poisson's ratio (ν) from the wave velocities.
The friction angle (ϕ') and cohesion (c') have been derived directly from direct shear tests. The material properties are reported in Table 4. Furthermore, in order to simulate the soil's damping characteristics in cyclic loading, Rayleigh damping has been defined. The damping ratio ξ has been taken to equal 10% [42]. The values of Rayleigh coefficients α and β are 2.04552 and 0.00434, respectively.

Structural Elements and Loads
In the finite element analysis, structural elements have been simulated by means of the elastic constitutive model.
The piles have been modelled as embedded beams. The embedded beam approach was introduced by Sadek and Sharhour [43]. In PLAXIS 3D, the embedded beam element is modelled as a three-node line element (beam element) that can cross a 10-node tetrahedral element, representing the soil in any place with any arbitrary orientation [44]. Due to the presence of the beam element, three extra nodes are introduced inside the 10-node tetrahedral element [45] (Figure 20). s 2020, 10, x FOR PEER REVIEW 17 of 32 The friction angle (φ') and cohesion (c') have been derived directly from direct shear tests. The material properties are reported in Table 4. Furthermore, in order to simulate the soil's damping characteristics in cyclic loading, Rayleigh damping has been defined. The damping ratio ξ has been taken to equal 10% [42]. The values of Rayleigh coefficients α and β are 2.04552 and 0.00434, respectively.

Structural Elements and Loads
In the finite element analysis, structural elements have been simulated by means of the elastic constitutive model.
The piles have been modelled as embedded beams. The embedded beam approach was introduced by Sadek and Sharhour [43]. In PLAXIS 3D, the embedded beam element is modelled as a three-node line element (beam element) that can cross a 10-node tetrahedral element, representing the soil in any place with any arbitrary orientation [44]. Due to the presence of the beam element, three extra nodes are introduced inside the 10node tetrahedral element [45] (Figure 20).  An elastoplastic model is used to describe the behavior of the interface. The interaction could involve a skin resistance as well as a foot resistance. For the skin and the base resistance, a failure criterion is used to distinguish between elastic and plastic interface behavior. For the interface, to remain elastic, the shear force (t s ) at a particular point is given by: |t s | < T max (6) For the plastic behavior: |t s | < T max (7) where T max is the equivalent local skin resistance at that point. The total pile bearing capacity, N pile , is given by: where F max is the base resistance and F skin is the axial skin resistance [46,47]. The base resistance has been calculated using the following equation.
where A p = base area of pile, c u = undrained cohesion at the base of the pile, N c = bearing capacity factor, assumed to equal 9, and σ v0,p = vertical total stress at the base of the pile. The skin resistance has been related to the strength of soil in which the pile is located: c i = R inter c soil (11) tan ϕ i = R inter tan ϕ soil (12) where τ i is the local shear stress resistance of the interface, ϕ i and c i are the friction angle and the cohesion of the interface, ϕ soil and c soil are the friction and the cohesion of the correspondent soil layer, R inter is the strength reduction factor, and σ' n = (σ' 2 + σ' 3 )/2 is the normal stress. The interaction has been modeled by choosing a suitable value for the strength reduction factor in the interface. R inter has been assumed to be 2/3.
For elasticity, the properties required are the Young modulus E and the unit weight of the beam element γ. Regarding geometric properties, the value of the cross-section area A and the moments of inertia have been calculated from the diameter of the massive circular pile. The embedded beam diameter determines the size in which plastic soil behavior is excluded.
Material properties of piles are reported in Table 5. Finally, the displacement and rotation at the pile head are both coupled with the displacement and rotation of the structural element in which the pile top is located.
The octagonal foundation plate has been modeled as Plate. After meshing, plates are composed of six-node triangular plate elements with six degrees of freedom per node. A plate has the following general and stiffness properties: thickness d, unit weight γ, Young's modulus in the first axial direction E 1 , Young's modulus in the second axial direction E 2 , in-plane shear modulus G 12 , out-of-plane shear modulus related to shear deformation over the first direction G 13 , out-of-plane shear modulus related to shear deformation over the second direction G 23 and Poisson's ratio ν 12 . The material behavior has been considered as orthotropic. Therefore: Table 6 shows the material properties of the plate. Finally, the octagonal pedestal has been modelled with an elastic constitutive model (E = 32,587,468 kN/m 2 , ν = 0.2, γ = 25kN/m 3 ).
In this study, operating and earthquake loads have been applied to the model. Operating dead load, D o , is the empty weight of equipment plus the maximum weight of contents (fluid load) during a normal operation [5]. It has been simulated as a cylindrical column with a total weight of 294 kN.
With regard to the earthquake loads, the seismic event that occurred in January 1693 has been chosen as a scenario earthquake. It was evaluated assuming the source to be along the Hyblean-Maltese fault and generating the 1693 seismic ground motion scenario [48][49][50]. The seismogram has been scaled to the value of 0.350 g, which corresponds to a return period of 975 years in the actual Italian seismic code NTC 2018 [51] (Figure 21). The predominant frequency of the 1693 input motion is 3.9 Hz. Furthermore, the response analysis to the earthquake on 13 December 1990 has been obtained using as base excitation as the record at the Sortino E-W component ( Figure 22) being the most representative of the earthquake. It has a PGA of 0.1 g and a predominant frequency of 1.82 Hz. Finally, in order to perform a parametric study with different input motions, the response analysis was also carried out using the recorded accelerogram of the 2012 earthquake at the Mirandola (MRN) station as excitation at the base of the model. On 20 May 2012, an earthquake of magnitude M W = 6.1 struck the Emilia Romagna Region of Italy and a small portion of the Lombardia region caused 27 deaths and considerable damage [52]. The MRN station is located at an epicentral distance of 13 km and is classified as a C site [42]. Therefore, a deconvolution analysis has been performed to obtain the appropriate input motion at the bedrock (PGA = 0.30 g, f p = 1.93 Hz) (Figure 23).   The schematic geometric model used to investigate the seismic response of the fullycoupled system is reported in Figure 24.   The schematic geometric model used to investigate the seismic response of the fullycoupled system is reported in Figure 24.   The schematic geometric model used to investigate the seismic response of the fullycoupled system is reported in Figure 24. The schematic geometric model used to investigate the seismic response of the fullycoupled system is reported in Figure 24.

Results
Regarding the two different alignments reported in Figure 19, the dynamic response has been investigated including its accelerations, response spectra, and amplification functions.
The numerical results from the analyses are shown in Figure 25 describing maximum accelerations with a depth using 1990, 2012, and 1693 seismic inputs.

Results
Regarding the two different alignments reported in Figure 19, the dynamic response has been investigated including its accelerations, response spectra, and amplification functions.
The numerical results from the analyses are shown in Figure 25 describing maximum accelerations with a depth using 1990, 2012, and 1693 seismic inputs. The values of the surface maximum accelerations and soil amplification factors R for the SSI and FF alignments are reported in Table 7. The results show that the 2012 and 1693 input motions are subjected to an evident amplification in the last 20 m and the 1990 accelerogram is subjected to an amplification in the last 10 m. The comparison between the alignment along the structure (SSI) and the alignment far from the structure (FF) shows that the presence of the structure generates a lower amplification. The soil amplification factors for FF alignment are greater than that achieved for SSI alignment.
In this case, the soil-structure interaction has beneficial effects. However, it is important to investigate and understand the SSI effects because the interaction can sometimes be detrimental. Furthermore, in both alignments, the soil amplification factors obtained using the 1990 and 1693 input motions are greater than the amplification value provided by Italian technical code [51], equal to 1.21 for soil type C. The smallest values of R were found for the 2012 accelerogram.
In Figure 26, the results are presented in term of response spectra at the surface, obtained by setting a structural damping of 5% and using the 1693 seismogram as inputs and the 1990 and 2012 accelerograms. For comparison, the elastic response spectrum provided by the Italian seismic code [51] is shown.
Considering the results obtained using the 1693 seismogram, for the free-field condition (Figure 26a), the maximum spectral acceleration Se,max = 1.94 g at T = 0.25 s is obtained, while, for the SSI alignment (Figure 26b), the maximum spectral acceleration The values of the surface maximum accelerations and soil amplification factors R for the SSI and FF alignments are reported in Table 7. The results show that the 2012 and 1693 input motions are subjected to an evident amplification in the last 20 m and the 1990 accelerogram is subjected to an amplification in the last 10 m. The comparison between the alignment along the structure (SSI) and the alignment far from the structure (FF) shows that the presence of the structure generates a lower amplification. The soil amplification factors for FF alignment are greater than that achieved for SSI alignment.
In this case, the soil-structure interaction has beneficial effects. However, it is important to investigate and understand the SSI effects because the interaction can sometimes be detrimental. Furthermore, in both alignments, the soil amplification factors obtained using the 1990 and 1693 input motions are greater than the amplification value provided by Italian technical code [51], equal to 1.21 for soil type C. The smallest values of R were found for the 2012 accelerogram.
In Figure 26, the results are presented in term of response spectra at the surface, obtained by setting a structural damping of 5% and using the 1693 seismogram as inputs and the 1990 and 2012 accelerograms. For comparison, the elastic response spectrum provided by the Italian seismic code [51] is shown.
In both cases, the maximum spectral acceleration is lower, considering full-coupled analysis. Furthermore, for a period T greater than 0.3 s, the elastic response spectrum provided by NTC 2018 [51] is more conservative than those obtained from FEM analysis.
Finally, considering the 1990 accelerogram, for the free-field condition (Figure 26a), the maximum spectral acceleration Se,max = 0.41 g at T = 0.23 is obtained, while, for the SSI alignment, the maximum spectral acceleration Se,max = 0.32 g at T = 0.25 is found.  Considering the results obtained using the 1693 seismogram, for the free-field condition (Figure 26a), the maximum spectral acceleration S e,max = 1.94 g at T = 0.25 s is obtained, while, for the SSI alignment (Figure 26b), the maximum spectral acceleration S e,max = 1.21 g at T = 0.24 s is found. A second less important period T = 0.58 s can be observed in the SSI condition.
The results obtained using the 2012 accelerogram show that the maximum spectral acceleration S e,max = 1.46 g at T = 0.25 s is obtained for the free-field condition (Figure 26a), while, for the SSI alignment (Figure 26b), the maximum spectral acceleration S e,max = 1.20 g at T = 0.25 s is found. A second less important period T = 0.05 s can be observed in FF and SSI conditions.
In both cases, the maximum spectral acceleration is lower, considering full-coupled analysis. Furthermore, for a period T greater than 0.3 s, the elastic response spectrum provided by NTC 2018 [51] is more conservative than those obtained from FEM analysis.
Finally, considering the 1990 accelerogram, for the free-field condition (Figure 26a), the maximum spectral acceleration S e,max = 0.41 g at T = 0.23 is obtained, while, for the SSI alignment, the maximum spectral acceleration S e,max = 0.32 g at T = 0.25 is found. Figure 27 shows the amplification functions A(f) evaluated as the ratio between the Fourier spectrum at the surface level and the Fourier spectrum at the bedrock level instead of the SSI.
Geosciences 2020, 10, x FOR PEER REVIEW 24 of 32 Figure 27 shows the amplification functions A(f) evaluated as the ratio between the Fourier spectrum at the surface level and the Fourier spectrum at the bedrock level instead of the SSI.
It is possible to observe that, between the amplification functions achieved for SSI and FF alignments, there are no substantial differences. Moreover, Figure 27 shows that the main resulting frequencies for SSI and FF conditions are equal to: fFF(I) = fSSI (I) = 2.1 Hz and fFF(II) = fSSI (II) = 3.7 Hz. Therefore, A(f) peaks move toward greater frequencies in comparison with the fundamental natural frequency of the soil (f0 = 1.37 Hz). These results are due to the higher predominant frequencies of the input motions.   It is possible to observe that, between the amplification functions achieved for SSI and FF alignments, there are no substantial differences. Moreover, Figure 27 shows that the main resulting frequencies for SSI and FF conditions are equal to: f FF (I) = f SSI (I) = 2.1 Hz and f FF (II) = f SSI (II) = 3.7 Hz. Therefore, A(f) peaks move toward greater frequencies in comparison with the fundamental natural frequency of the soil (f 0 = 1.37 Hz). These results are due to the higher predominant frequencies of the input motions.
The fundamental period of the fixed-base structure is equal to 0.54 s. Based on the response spectra obtained using the 1693 seismic input (Figure 26), a spectral acceleration of 0.99 g has been obtained, according to NTC 2018 [51], while lower values of 0.66 g (FF alignment) and 0.85 g (SSI alignment) have been found from FEM numerical analysis. Therefore, for the period of the structure under consideration, the NTC 2018 [51] spectrum is more conservative. Moreover, the soil-structure interaction leads to a designed acceleration greater than that required in the FF condition due to the second period T = 0.58 s in the SSI condition.
The period of the structure resting on the soil can be calculated as the ratio between the Fourier amplitude spectra computed at the top and at the bottom of the structure [41]. The resulting period of the structure, TSTRU, SSI, is equal to 0.95 s (fSTRU,SSI = 1.05 Hz). In this case, a spectral acceleration of 0.65 g has been obtained, according to NTC 2018 [51], while lower values of 0.37 g (FF alignment) and 0.40 g (SSI alignment) have been found from FEM numerical analysis. SSI has beneficial effects because the spectral accelerations are lower than those required for the period of the fixed-base structure. Figure 29 shows the distribution of the bending moment along the plate at the end of the dynamic phase using the 1693 seismic input. The results of the moment distribution indicate that the maximum bending moments are distributed based on the pile heads. The fundamental period of the fixed-base structure is equal to 0.54 s. Based on the response spectra obtained using the 1693 seismic input (Figure 26), a spectral acceleration of 0.99 g has been obtained, according to NTC 2018 [51], while lower values of 0.66 g (FF alignment) and 0.85 g (SSI alignment) have been found from FEM numerical analysis. Therefore, for the period of the structure under consideration, the NTC 2018 [51] spectrum is more conservative. Moreover, the soil-structure interaction leads to a designed acceleration greater than that required in the FF condition due to the second period T = 0.58 s in the SSI condition.
The period of the structure resting on the soil can be calculated as the ratio between the Fourier amplitude spectra computed at the top and at the bottom of the structure [41]. The resulting period of the structure, T STRU, SSI , is equal to 0.95 s (f STRU,SSI = 1.05 Hz). In this case, a spectral acceleration of 0.65 g has been obtained, according to NTC 2018 [51], while lower values of 0.37 g (FF alignment) and 0.40 g (SSI alignment) have been found from FEM numerical analysis. SSI has beneficial effects because the spectral accelerations are lower than those required for the period of the fixed-base structure. Figure 29 shows the distribution of the bending moment along the plate at the end of the dynamic phase using the 1693 seismic input. The results of the moment distribution indicate that the maximum bending moments are distributed based on the pile heads.
The seismic response of the piled foundation has been further examined in terms of the dynamic pile bending moment. In fact, the seismic response of piles involves the inertial interaction between the structure and the piles as well as the kinematic interaction between the piles and the soil foundation. Thus, for pile foundations built in seismic areas, the demands to sustain load and deformation during an earthquake will likely be the most severe in their design life.
FEM numerical analysis. SSI has beneficial effects because the spectral accele lower than those required for the period of the fixed-base structure. Figure 29 shows the distribution of the bending moment along the plate a the dynamic phase using the 1693 seismic input. The results of the moment d indicate that the maximum bending moments are distributed based on the pil During an earthquake there are two sources of loading on piles: (i) "inertial" loading of the pile head caused by the lateral forces imposed on the superstructure, and (ii) "kinematic" loading along the length of the pile caused by the lateral soil movements developed during the earthquake, assuming zero inertia at the superstructure [12,61,62] As a consequence, piles should be designed for the following two loading conditions: (a) inertia forces on the superstructure transmitted on the heads of the piles in the form of axial and horizontal forces and moment; (b) soil deformations arising from the passage of seismic waves, which impose curvatures and, thereby, a lateral strain on the piles along their whole length.
Kinematic bending moments are significant especially with regard to stiff pile caps and soil layer interfaces characterized by a contrast of stiffness.
While there is ample experience for carrying out the equivalent static analysis for inertial loading (type (a)), no specific method or procedures are available to predict deformations and the bending moment for kinematic loading (type (b)).
With this aim, the computed bending moments distribution along the pile length are reported in Figure 30, for pile 4, as an example. In this case, the numerical simulation provides the bending moment distribution along the pile length related to inertial loading and kinematic lateral soil movements. For the 1693 seismic input, a bending moment of 134 kNm can be observed at the pile head, characterized by a fixed connection between the plate and the piles head.
The maximum bending moment (365 kNm) occurs at a depth of about 12.5 m, reducing progressively to zero close to the piles tip. Similar trends are obtained using the 1990 and 2012 seismic inputs ( Figure 30). These maximum bending moments can result in a kinematic interaction effect between the piles and the soil.
Kinematic pile bending tends to be amplified close to the interface ( Figure 31) between the softer alluvial deposits and stiffer sandy clay (Table 3), which typically occurs due to a kinematic interaction effect [12].    Kinematic bending of a free-head pile in a two-layer profile (After Castelli et al. [12], modified). Figure 31. Kinematic bending of a free-head pile in a two-layer profile (After Castelli et al. [12], modified).
Finally, the stress state in the finite element model along a vertical cross section of the mesh at the end of the dynamic phase (for the 1693 seismic input) is displayed in Figure 32 in terms of mean total stress p = (σ 1 + σ 2 + σ 3 )/3, where σ 1 is the largest compressive principal stress and σ 3 is the smallest compressive principal stress (σ 1 ≤ σ 2 ≤ σ 3 ). The variation of the stress state in the area close to the pile tips is likely related to a more pronounced interaction effect of the piles in the group with respect to the fixed connection of the top.

Conclusions
In the city of Augusta, located in Southeastern Sicily (Italy) and characterized by a high seismic hazard, there are important, industrial, petrochemical facilities. They are classified with high risk according to the Seveso III Directive based on major-accident hazards involving dangerous substances [30].
The seismic vulnerability of petrochemical facilities is often analyzed by applying codes and guidelines developed primarily in the design of new and existing buildings. On the contrary, the evaluation of seismic safety of these facilities requires specific evaluation issues.
The idea of the proposed study is to provide a practical guidance to people involved in the seismic design of petrochemical facilities, as a best practice for the vulnerability evaluation. It is intended to serve the following objectives: • to better understand the intent behind certain provisions of seismic design codes, so that they can be more properly and uniformly applied to structures and systems typically found in petrochemical facilities; • to provide background information on technical areas that are related to the seismic evaluation of petrochemical facilities; • to provide specific guidance to the seismic evaluation of petrochemical facilities; • to provide practical analytical guidance specifically applicable to petrochemical facilities.
With this aim, the dynamic response of a fully-coupled, soil-structure system has been carried out by means of a PLAXIS 3D numerical code.
For the site characterization of soil, deep site investigations have been performed. The response of the system has been analyzed, considering and not considering the Soil-Structure Interaction (SSI), which demonstrates that the SSI can lead to alterations of the free-field motion. In the literature, only a limited number of studies have been carried out considering deep foundations embedded in multi-layered soil profiles. This is likely due to the high computational cost and the skills required to carry out numerical simulation using Finite Element analysis. Therefore, in our opinion, a study like the one proposed in the paper should be welcomed. In fact, FEM is a very viable method to capture the real behaviour of a soil-pile interaction in the three-dimensional domain and a continuum-based analysis can be used to simulate, in a realistic way, the kinematic pile behaviour under a seismic excitation.

Conclusions
In the city of Augusta, located in Southeastern Sicily (Italy) and characterized by a high seismic hazard, there are important, industrial, petrochemical facilities. They are classified with high risk according to the Seveso III Directive based on major-accident hazards involving dangerous substances [30].
The seismic vulnerability of petrochemical facilities is often analyzed by applying codes and guidelines developed primarily in the design of new and existing buildings. On the contrary, the evaluation of seismic safety of these facilities requires specific evaluation issues.
The idea of the proposed study is to provide a practical guidance to people involved in the seismic design of petrochemical facilities, as a best practice for the vulnerability evaluation. It is intended to serve the following objectives: • to better understand the intent behind certain provisions of seismic design codes, so that they can be more properly and uniformly applied to structures and systems typically found in petrochemical facilities; • to provide background information on technical areas that are related to the seismic evaluation of petrochemical facilities; • to provide specific guidance to the seismic evaluation of petrochemical facilities; • to provide practical analytical guidance specifically applicable to petrochemical facilities.
With this aim, the dynamic response of a fully-coupled, soil-structure system has been carried out by means of a PLAXIS 3D numerical code.
For the site characterization of soil, deep site investigations have been performed. The response of the system has been analyzed, considering and not considering the Soil-Structure Interaction (SSI), which demonstrates that the SSI can lead to alterations of the free-field motion.
The results of the numerical modeling have also been compared with those suggested by the Italian Technical Regulations [51]. The findings of this study show that many factors have to be considered on the seismic design of new structures and/or apparatuses or retrofitting of existing ones. The main considerations can be summarized as follows.

•
The FEM analyses show a strong amplification in the last 20 m up to ground level for both Free-Field (FF) and soil-structure interaction (SSI) alignments using the 2012 and 1693 seismic inputs. The 1990 accelerogram induces a larger amplification in the upper 10 m of soil. Moreover, the presence of the structure generates a lower amplification as compared to the free-field condition; • The influence of the soil-structure interaction is also indicated by reducing the maximum spectral acceleration at the main period of T = 0.25 s. However, a second, less important, period T = 0.58 s, corresponding to a spectral acceleration of 1.0 g, can be observed in the SSI condition using the 1693 seismogram; • Taking into account the SSI effects, a beneficial effect can be observed. The spectral accelerations obtained when considering the structure resting on the soil are lower than what is required for fixed-base structure's period; • The main resulting frequencies for the SSI condition are equal to: f SSI (I) = 2.1 Hz and f SSI (II) = 3.7 Hz. They are far from the frequency of the structure f STRU,SSI = 1.05 Hz.

•
The soil amplification factors R derived for SSI and FF alignments using the 1990 and 1693 input motions are greater than the amplification value provided by the Italian Technical Code [51]. Instead, lower values are obtained using the 2012 accelerogram. The comparison between the elastic response spectra obtained by numerical analyses using the 1693 and 2012 seismic inputs and the same provided by NTC 2018 shows that FEM spectral accelerations are lower than those given by NTC 2018 for periods greater that 0.30 s. In particular, for the period of the structure under consideration (T STRU,SSI = 0.95 s), the spectral acceleration given by NTC 2018 is equal to 0.65 g.

•
To investigate the effect of the kinematic and inertial interaction on pile bending moments, the distribution of the peak pile bending moments has been studied. It can be observed that the maximum bending moments occur at a depth of about 12.5 m, and it can result in a kinematic interaction effect between the piles and the surrounding soil.
As a consequence of this last indication, it can be stated that, if kinematic effects are ignored and only inertial effects are considered, the maximum bending moment distribution along the piles length can be underestimated.