A risk‐targeted seismic performance assessment of elephant‐foot buckling in walls of unanchored steel storage tanks

Unanchored steel storage tanks are used in industrial plants to store oil and other petrochemical liquids. In a major earthquake, such tanks may sustain different types of damage, including elasto‐plastic buckling in the shape of an elephant's foot, which is the object of the research. The methodology introduced in this paper consists of cloud‐based 1D site response analysis (SRA) and pushover‐based seismic performance assessment of the tank wall using a 3D non‐linear model of the tank. The SRA, which considers the recorded ground motions on the rock outcrop and the sample of interval shear‐wave velocity profiles, is carried out to estimate spectral accelerations at the tank impulsive period. By adopting the code‐based model of hydrodynamic pressures, the seismic demand on the tank wall is then calculated by pushover analysis. The risk‐targeted decision model is used for safety verification. The proposed methodology is demonstrated through an example of a large liquid storage tank for which it is shown that the fragility function for peak ground acceleration at the rock outcrop is on the left side of the target fragility function, indicating that the tank cannot be considered safe from elephant‐foot buckling. The introduced assessment process is relatively simple if the cloud‐based SRA is automated. However, further experimental and numerical investigations are required to confirm the validity range of the pushover analysis for the seismic performance assessment of tanks.


NOVELTY
• A risk-targeted seismic performance assessment of elephant-foot buckling (EFB) of steel storage tanks is introduced.• The methodology decouples the seismic performance assessment of EFB to 1D site response analysis (SRA) and non-linear pushover analysis of the refined model of the tank.• The risk-targeted intensity measure (IM)-based decision model is incorporated into the introduced methodology.

INTRODUCTION
Above-ground liquid storage tanks are used in industrial facilities to store and process a variety of hazardous liquids and liquid-like materials including oil, liquefied natural gas, chemical fluids and waste liquids.Several major earthquakes around the world have demonstrated inadequate seismic safety in many liquid storage tanks.2][3] Post-earthquake field inspections have indicated distinct failure modes of the tanks.Damage to tank components was most commonly related to elephant-foot buckling (EFB) of the wall of an above-ground liquid storage tank (Figure 1), roof damage, base plate failure and nozzle failure. 45][6][7][8] The occurrence of NaTech events should thus be limited, which is often addressed by quantitative risk analysis aimed at preventing earthquake-induced leakage with high reliability to avoid environmental and human harm. 91][12][13][14] In the framework of quantitative risk analysis of liquid storage tanks, fragility functions can be defined in the form of probit functions. 15However, damage states may be defined by engineering demand parameters (EDPs) obtained from the seismic analysis.Few recent studies have considered the effect of base uplifting on the tank wall and addressed the seismic response of the tank. 7,16,17Malhotra and Veletsos 10 studied the seismic response of tanks by developing a lumped mass-spring model.They decoupled the impulsive and convective liquid motion, indicating a vast difference in their periods.Their model has a rigid base plate that is subjected to uplifting during seismic loading.They later added a rotational spring representing rocking resistance to the tank base using the beam method. 10Recently, the Malhotra and Veletsos model was used for risk assessment. 11Vathi and Karamanos 16 considered the non-linear effects of base uplifting on the tank wall and redefined the rotational spring mechanical behaviour based on non-linear static analysis of a detailed finite element (FE) model.Recent studies have focussed on better defining the non-linear characteristics of the base plate using simplified 3D models. 7,17In these models, the nonlinear behaviour of the tank is concentrated in springs at the edge of the tank connecting the base plate to the ground.Few recent experimental studies have focussed on the stress distribution on the tank wall in above-ground liquid storage tanks. 18,19any researchers 12,13,20 studied above-ground liquid storage tanks using refined FE models.Their studies focussed on two models applied to the liquid domain: an arbitrary Lagrangian-Eulerian model 12 and an acoustic model. 13,20Although wall buckling may lead to leakage and loss of tank containment, non-linear effects or large deformations affecting the tank wall due to EFB have not been thoroughly studied. 4,6However, it was observed that the occurrence of EFB does not necessarily trigger the loss of containment. 5,6,8 Manos and Clough (1983) presented tanks where EFB occurred, but the leakage was not observed in some cases.In their post-earthquake observations, the EFB-based leakage occurred at the pipe connections and bolted joints. 5Therefore, the presence of discontinuities on the tank wall (e.g., piping, nozzles, manholes) and large strains in the area where EFB takes place increase the probability of rupture (and thus loss of content).In this respect, computationally efficient tools capable of capturing the non-linear behaviour of the tank wall are needed.0][21] Such an analysis method may become computationally too demanding if used in a probabilistic approach. 22Thus, developing simplified procedures to assess the seismic performance of tanks considering F I G U R E 1 Elephant-foot buckling of an above-ground storage tank. 60e non-linear behaviour of the tank wall is sensible.Substantial efforts have been made in developing simplified models of above-ground liquid storage tanks.19,21 The seismic performance of structures depends on the ground motion intensity, as was also observed for storage tanks.7,12,16,20 It is known that ground motion intensity at the site is, to some extent, affected by the site properties 23 Passeri et al. 24 reported that the soil deposit properties affect the site response and the seismic action at the ground surface and thus, for example, at the structure foundation elevation.However, the authors used different ground motion intensity measures (IMs) for seismic hazard study.Peak ground acceleration (PGA) at reference rock is often used 7,8,25,26 as it is always available from the probabilistic seismic hazard analysis that is required in conventional risk analysis (e.g., [27][28][29][30] ).Some other authors used the spectral acceleration at the tank's impulsive period, which is more directly related to the seismic response of the tank.Consequently, it is an efficient ground motion IM for seismic fragility analysis.8,[31][32][33] However, the seismic hazard function for spectral acceleration at the tank's impulsive period and at the elevation of the tank's foundation is often unavailable. In th case, the site response analysis (SRA) can be performed to convert ground motion IMs used in seismic response analysis and seismic hazard/risk analysis.
The objective of the research presented in this paper is the simplified risk-targeted seismic performance assessment of the EFB limit state in the wall of steel broad storage tanks.The proposed methodology involves 1D SRA, pushover-based seismic performance assessment of the tank wall using a refined FE model of the tank and the risk-targeted IM-based decision model.The 1D SRA provides the acceleration response spectra at the tank foundation and, consequently, the relationship between the ground motion IM considered in seismic hazard analysis and that considered in the seismic response analysis.The pushover analysis is used to couple the seismic demand related to the spectral acceleration at the impulsive period S e with the limit-state capacity defined by the stresses of the tank wall.
The remainder of the paper is organised as follows.The proposed methodology is described in Section 2. A seismic performance assessment of an existing unanchored steel storage tank is presented in Section 3. The site properties at the tank location, the ground motions and the results of the SRA are presented in Section 3.1.The tank properties, a 3D nonlinear model, the loading on the tank, the pushover analysis and assessment of the seismic performance of the tank walls are presented in Section 3.2.The limitations of the methodology are discussed in Section 4. The findings of the study are addressed in Section 5.

DECOUPLED METHODOLOGY FOR SEISMIC PERFORMANCE ASSESSMENT OF ELEPHANT-FOOT BUCKLING OF WALLS OF UNANCHORED STEEL STORAGE TANKS
The proposed methodology decouples the problem into two independent processes: (A) SRA and (B) pushover-based seismic performance assessment of the EFB of the tank wall, as shown in Figure 2. The SRA consists of definition of the site properties (A1), ground motion selection (A2) and equivalent linear 1D site response analyses (A3) 34 conducted in the DEEPSOIL open-source platform. 35The effects of the intrinsic spatial variability of soil deposit properties are systematically considered by generating a sample of shear-wave V s soil profiles based on the stochastic modelling of 1D soil profiles recently proposed by Passeri et al. 24 Step A1 in the proposed methodology results in a sample of the V s soil profiles, modulus reduction and damping curves.The 1D SRA solves the problem of vertical propagation of horizontal shear waves of an earthquake from the bedrock level through a horizontally layered site.In this regard, only the ground motions recorded on the rock outcrop are considered.
7][38] Site response analyses (A3) are performed for all 1D soil models defined by the sample of V s soil profiles (A1) and for all selected ground motions (A2).

F I G U R E 2
Flowchart of Processes A and B in the proposed methodology for elephant-foot buckling limit-state safety verification of unanchored liquid storage tanks.
Step A3 results in a sample of acceleration response spectra at the foundation level of the tank.The spectral acceleration at the impulsive period of the tank S e is an input parameter for estimating dynamic pressure and stresses and deformation on the tank wall, linking with Process B in the proposed methodology.
Process B includes pushover analysis based on the refined FE model of the tank, developed from non-linear shell elements capable of simulating material non-linearity and supporting large strains (B1).Adequate interface elements, allowing for uplifting of the tank base plate from the rigid foundation, are also used.Thus, the geometric non-linearity at the tank-foundation interface due to base uplifting is considered.The dead loads on the tank include the gravity load of the structure and the hydrostatic pressure (B2).The seismic action on the tank is introduced by a force-controlled pushover analysis (B3).The load pattern on the tank used in pushover analysis is based on the hydrodynamic pressure models due to the horizontal and vertical components of the ground motion 32 (see Section 2.2).As the shape of the hydrodynamic pressure is considered constant for a given tank configuration, only one load pattern is defined for the pushover analysis.
The pushover analysis of the tank is performed using FEM software such as Abaqus 39 via the Riks method, which is suitable for analysing stability problems.The Riks method simultaneously solves the non-linear static problem for loads and displacements along the static equilibrium path in load-displacement space. 40The pushover analysis provides the relationship between the base shear and the maximum axial compressive stress acting on the tank wall σ, which is used as the EDP to define the occurrence of EFB.The pushover analysis is performed until non-linear EFB effects.The base shear from pushover analysis is equal to the resultant of the hydrodynamic pressure imposed on the tank wall.As the code-based hydrodynamic pressure model of the tank depends on the spectral acceleration at the impulsive period S e , the base shear from the pushover analysis can be expressed in terms of S e .
The code-based hydrodynamic pressure model allows coupling of the SRA, which provides the demand in terms of S e , and the pushover analysis to obtain the seismic demand for EFB in terms of the maximum axial compressive stress in the tank walls (Section 2.2).Coupling between the SRA (A3) and pushover analysis (B3) is performed in step B4 of the proposed methodology.However, S e is converted to the corresponding PGA at the rock outcrop because is considered the seismic IM often used to study the seismic vulnerability of tanks. 7,8,25,26The PGA at the rock outcrop is the ground-motion IM for development of the seismic demand model (B4).The preliminary result of step B4 is a cloud of seismic demand points expressed in terms of PGA at the rock, and σ, the EDP.These pairs represent an input for development of the seismic demand model, which for cloud analysis is based on linear regression analysis.

Site response analysis (Process A)
SRA is conducted in three independent steps.The definition of site properties (A1) begins with defining the reference soil profile of shear-wave velocity V s and ends with the sample of the interval V s profiles paired with modulus degradation and damping curves.SRA involves the model for the shear-wave velocity profile.The model can be based on a reference site profile obtained from in situ soil measurements at a specific location. 24In the absence of detailed soil measurements, the site categorisation in Eurocode 8 41 can be used, which is based on the depth of the soil layer corresponding to the seismic bedrock formation H 800 and the harmonic average velocity in the top 30 m of soil layers V sS, 30 .In this case, the reference shear-wave velocity profile can be based on the empirical formula for the maximum shear modulus G max . 42,43 max () = 1000 2,max ( ′  ()) 0.5 where σ ′ m is the mean effective confining stress dependent on the depth z (σ′ m = ρ g z), and K 2,max is a constant dependent on the soil density ρ. 42 The reference V s profile of the site can be calculated based on the theory of 1D shear-wave propagation: The soil profile in Equation ( 2) is a continuous profile of shear-wave velocities V S (zZ) discretised into n layers, an interval reference V s profile used in the 1D SRA.The effects of the intrinsic spatial variability of soil deposit properties are systematically considered according to the stochastic modelling recently proposed by Passeri et al. 24 This is achieved by generating a sample of interval V s profiles, referred to as equivalent V s profiles.
The randomisation approach proposed by Passeri et al. 24 is based on the interval reference V s .It starts by generating the layer interfaces (the depths at the bottom of the soil layers measured from the ground surface) by assuming a nonhomogenous Poisson process 44 that considers the depth-dependent rate of layer interfaces (number of layer interfaces per meter) λ(z) =  1 (+ 2 )  3 , with constants c 1 , c 2 and c 3 proposed by Toro. 44In the first step, the n layer thicknesses are generated by the inverse method, simply assuming a homogenous Poisson process: where d is the layer thickness based on a unit exponential distribution (λ = 1), and F(d) is the cumulative distribution function.The n layer thicknesses generated according to Equation (3) are transformed in the second step into actual layer thicknesses considering the depth-dependent rate of layer interfaces.This transformation requires the inverse of the cumulative rate function based on λ(z), previously derived 44,45 where u is the sum of the randomly generated layer thicknesses d up to a given layer.Λ −1 (u) represents the layer interfaces for a given layer, used to calculate the layer thicknesses t i .
The positions of the interfaces in the randomly generated profile should be bounded by the characteristics of the interval reference V s profile.The boundary criteria for sampling the layer thickness are defined by the minimum and maximum layer thickness, which are considered to be equal to the minimum and maximum layer thickness of the interval reference V s profile, respectively.If the interval reference V s profile is not obtained from in situ soil measurements at a specific location, it is recommended that the interval reference V s profile layers have a minimum thickness of 1 m and a maximum thickness of 20 m.Further, it is recommended that the generated layer thickness is bounded by an upper multiplication factor of 1.25 and a lower multiplication factor of 0.75 with respect to the layer thickness of the interval reference V s profile. 24he soil-column randomisation approach generates a sample of interval V s profiles using the travelling time of the shear wave from the soil depth z to the ground surface tt S,Z .Assuming tt S,Z is a log-normally distributed variable, 24 a random profile of tt S,Z is generated by solving the left side of Equation ( 5) with respect to tt S,Z , where the mean value  , is obtained from the interval reference V s profile considering Equation (6).Equation (6) shows the correlation between the interval shear-wave velocity V s , the harmonic average shear-wave velocity V S,Z , and the cumulated travel time tt S,Z , considered as  , in the reference interval V s profile.A random tt S,Z profile is calculated at the mid-elevation of each layer.Thus, in Equations ( 5) and ( 6), the variables with index i or dependent on i refer to the ith layer.The logarithmic standard deviation  ( , ) in Equation ( 5) can be assumed to be 0.015. 24The right side of Equation ( 5) is used to compute the S(i) for ith layer.In this context, ρ(i) is the interlayer correlation coefficient, S(i-1) is a normal random variable for i-1th layer, and ε is a random variable with zero mean and unit standard deviation.The interlayer correlation coefficient ρ(i) was proposed by Toro 44 and is a function of z corresponding to the mid-elevation of the ith layer and distance between layer mid-points r. 24 The variables z and r are defined using the randomly generated profile obtained from Equation (4).
Once the tt S,Z profile is generated by Equation ( 5), it is converted to a V S,Z profile in the left side of Equation ( 6).However, the V S,Z profile does not have a direct application in the 1D SRA.Thus, it is converted to the corresponding interval V s profile.The shear-wave velocity of the ith layer is calculated using Equation ( 6) based on the known depth, layer thicknesses and V S,Z .The proposed randomisation process of the interval V s profile avoids double-counting of uncertainties that could occur if the soil profile is randomised directly to the interval V s profile (space-and time-dependent) and the depth of the interfaces (spatial variable).The randomisation approach is repeated as many times as necessary.
The Passeri et al. 24 model accounts for the variability of the bedrock shear-wave velocity.As with the soil deposit, randomisation of the bedrock shear-wave velocity  ℎ  is executed by means of a standard normal random variable  ℎ  : where  ℎ is the Pearson's correlation coefficient indicating the correlation between the bedrock depth  ℎ and the shearwave velocity at bedrock  ℎ  ; its value is considered to be 0.508. 24 ℎ is a standard normal random variable for bedrock depth  ℎ (Equation 7), and ε is a random variable with zero mean and unit standard deviation. ℎ can be calculated using Equation (8), as the randomised bedrock depth  ℎ derives from the soil deposit layering randomisation in a physically based approach. ℎ and  ℎ  are the depth and shear velocity, respectively, characterising the bedrock in the interval reference V s profile;  ln( ℎ ) is the standard deviation of natural logarithms of bedrock depth assumed to be 0.090. 18The randomised bedrock shear-wave velocity is evaluated using Equation ( 7) assuming the logarithmic standard deviation  ln( ℎ  ) is 0.101. 24n addition to the sample of interval V s profiles, the modulus degradation and damping curves are defined to complete the soil model (A1).In this study, the modulus degradation and damping curves were taken from the literature. 46 Darendeli 46 reviewed soil tests conducted at the University of Texas at Austin and proposed a four-parameter soil model that describes the change in normalised shear modulus G and material damping ratio D 46 : Masing +  min (10)   where γ r and a are the reference strain and curvature coefficient, respectively, dependent on the soil type and loading conditions; G max is the shear modulus at zero shear strain γ from the generated equivalent V s profile; D min is the small-strain material damping ratio; b is the scaling coefficient and D Masing (a, γ r , γ) is the damping estimated from the Masing behaviour, dependent on a, γ r and γ.
The modulus degradation and damping curves are defined by the hyperbolic model with non-Masing behaviour (MRDF) included in DEEPSOIL. 35This model applies the shear modulus degradation and damping ratio curves as functions of shear strain. 47The shear modulus degradation and damping ratio curves from Equations ( 9) and ( 10) 46 are simultaneously fit by introducing a reduction factor that alters the Masing rules. 48n the next step of Process A, the ground motions at the rock outcrop for the 1D SRA are defined (A2).As the SRA is based on the cloud analysis, the ground motions were selected from strong ground motion databases considering the magnitude, distance and V s,30 intervals. 49Such an approach implies no need for ground motion scaling.However, due to the lack of strong ground motions recorded on a stiff or rock site, the selected ground motion can be further scaled to make the ground motion set representative of the entire range of ground motion IMs.The final step of Process A is execution of the site response analyses (step A3) in DEEPSOIL, 35 which supports non-linear and equivalent linear analysis.In this study, equivalent linear analysis in a frequency domain was used.It considers the non-linear hysteretic response through the Kelvin-Voigt model, 34,35,42 which is iteratively modified at each soil layer with consideration of the modulus reduction and damping curves and the effective shear strain at a given soil layer.An effective shear strain is computed for each soil layer after the analysis considering a given ground motion and an initial estimate of shear modulus and critical damping ratio.The effective shear strain usually corresponds to 65% of the peak strain. 34odulus degradation and damping curves are then used to update the shear modulus and damping coefficient in each soil layer.An iterative scheme is used to produce a converged solution. 50o guarantee the prediction accuracy, in the equivalent linear 1D SRA, the magnitude of the shear strain γ, the parameter most influencing the analysis accuracy, 51 is usually limited to γ max < 0.35%.The number of analyses is equal to the number of generated equivalent V s profiles obtained in step A1 multiplied by the number of ground motions from step A2.The cloud-based SRA results in the acceleration response spectra for each soil profile and each ground motion.Based on the acceleration response spectra, the sample of the spectral acceleration at the impulsive rigid period of vibration of the tank-foundation system S e is considered in defining the hydrodynamic pressure on the tank 32 and the EDP using pushover analysis, explained as follows.The relationship between S e and the PGA values at the rock outcrop from input ground motions is considered because the PGA at the rock outcrop is the seismic IM considered in the probabilistic seismic hazard analysis and the resulting seismic hazard curve.PGA is often used in developing seismic demand models of liquid storage tanks. 25,26

Pushover-based seismic performance assessment (Process B)
The refined FE model of the tank is used for pushover analysis.In step B1, the tank properties and 3D model are defined.In this study, the tank configuration for pushover analysis is limited to broad unanchored steel storage tanks and considers the effects of uplifting and wall buckling phenomena on seismic vulnerability.As only broad tanks are considered, they can be classified according to the draft of the new Eurocode 32 as rigid tanks concerning the impulsive component of their response.The criteria for the ratio of the filling height H and tank radius R, and the ratio of the tank radius R and bottom wall thickness t are presented in table A.10 in the draft of the new Eurocode. 32Along the tank wall, the wall thickness decreases from bottom to top due to a reduction in the liquid pressure acting on the wall. 52n the scope of this study, the tank is considered filled to the maximum capacity level.Thus, the seismic performance of broad unanchored steel storage tanks with a filling height H = 0.9H tot was taken into account. 7,53 I G U R E 3 (A) gravity load (W) and hydrostatic pressure (p s ) applied through non-linear static analyses and (B) hydrodynamic pressures on the tank due to horizontal seismic action  ,ℎ and due to vertical seismic action  , , applied through pushover analysis.
The 3D FE model of the tank (Abaqus 39 ) uses non-linear shell elements capable to catch large deformations, and the non-linear stress-strain relationship of the tank material (B1).The 3D non-linear modelling also involves the nonlinearities due to the base uplifting and sliding (see Section 3).For simplicity, the foundation is considered to be rigid, and the soil-structure interaction is considered in the estimation of the period of the impulsive mode as defined in the draft of the new Eurocode. 32he loads on the tank are defined in step B2.The dead load is specified as the gravity load of the structure and the hydrostatic pressure p s .The load pattern imposed in pushover analysis is based on the hydrodynamic pressure models defined in Equations ( 11) and (12). 32The load pattern for pushover analysis is based on the patterns of the hydrodynamic pressure due to horizontal seismic action p ir,h and vertical seismic action p ir,v , 7,8,54-58 which are applied on the tank accounting for the SRSS rule, 32

√
,ℎ 2 +  , 2 .The p ir,h and p ir,v are defined as follows 32 : where   is the liquid density, H is the filling height,  = ∕ is the dimensionless height and z is the elevation measured from the bottom to the top of the tank. ,ℎ (, ) is the dimensionless impulsive rigid pressure function, dependent on the slenderness ratio of the tank ( = ∕) and dimensionless height .Γ ,ℎ is the participation factor of the impulsive rigid pressure component,  is the circumferential angle,  , is the spectral acceleration in the vertical direction and   is the spectral acceleration in the horizontal direction at the impulsive rigid period considering the soil-foundation-tank interaction  ,ℎ .For simplicity, the  , is considered the same as the spectral acceleration in the horizontal direction   (i.e.,  , =   ).Thus the stresses in the tank obtained from the pushover analysis account for the hydrodynamic pressure due to the horizontal and vertical ground motion components (Equations 11 and 12).However, the seismic IM for the fragility and risk analysis refers only to the horizontal ground motion component.The loads are applied to the tank model in three loading steps.In the first and second steps, the tank is loaded with the gravity load and the hydrostatic pressure of the liquid (Figure 3A), respectively, through non-linear static analyses.The hydrodynamic pressure is applied in the third step through incremental non-linear static (pushover) analysis with large displacements (Figure 3B).
The pushover analysis of the refined FE tank model (step B3) is performed by FEM software (e.g., Abaqus 39 ) at least until the formation of elasto-plastic EFB at the bottom of the tank's wall.The EFB limit state can generally be defined by a deformation criterion of the steel in the region of EFB.However, the EFB limit state is often defined by a stress criterion. 6,7,11,25,59A similar approach was proposed by Rotter 60 and is considered in the Eurocode 61 : where σ m is the maximum axial compressive stress corresponding to the EFB limit state and p ir is the maximum interior pressure acting on the tank resulting from the hydrostatic pressure p s and the hydrodynamic pressure , which accounts for horizontal and vertical ground motion components based on the SRSS rule. 32e σ cl is the ideal critical buckling stress dependent on tank wall geometry parameters R and t.The f y and E are the yield stress and elastic modulus of the tank wall material, respectively, and r = (R/t)/400.In the verification of the EFB limit state according to Equation ( 13), the axial compressive stress in the tank wall was estimated directly from the pushover analysis.
As the code-based hydrodynamic pressure models of the tank (Equations 11 and 12) depend on the spectral acceleration at the impulsive period S e , the S e -EDP relationship is obtained directly from pushover analysis.According to the Eurocode, 61 the axial compressive stress is considered an effective EDP to define the EFB limit state.Consequently, the spectral acceleration at the impulsive period S e and axial compressive stress σ is obtained from the pushover analysis, as demonstrated in step B3 in Figure 2.Alternatively, the code-based model for estimating the seismic demand for the axial compressive stresses in the tank's wall can be replaced by performing dynamic analyses based on a simplified model of the tank, 7,8,25 which directly account for the effect of the entire ground motion response history.
Generally, S e can be used as the IM in the seismic demand model.However, the seismic hazard analysis is often available only for the rock outcrop.In addition, the PGA is often considered in studies on seismic vulnerability of tanks. 25,26Thus, it was decided to link S e at the level of the tank foundation to the PGA at the rock.This can be done based on the input and output acceleration response spectra from the SRA.Thus, the intermediate result of step B4 is a cloud of seismic demand points (B4 in Figure 2) expressed in terms of PGA at the rock, and the axial compressive stress σ.These points represent an input for development of the seismic demand model, which in cloud analysis is based on linear regression analysis of the natural logarithm of IM and EDP. 62The probabilistic seismic demand model can be defined as 49 : where σ is the median demand in terms of the maximum axial compressive stress in the tank walls, a and b are regression coefficients, σ i is the maximum axial compressive stress in the tank wall for ith point from the cloud of points and n is the number of cloud points.Note that the choice of PGA as IM for the seismic fragility analysis is not a limitation of the proposed methodology.For example, the spectral acceleration at the tank's impulsive period at bedrock or the ground surface can be considered as IM for seismic fragility analysis.The intersection of the seismic demand defined in Equation ( 14) and the EFB limit state defined in Equation ( 13) defines the median PGA causing EFB based on the definition of the EFB limit state PGA EFB (B4 in Figure 2).However, given the uncertainty present in the definition of the EFB limit state, as discussed in ref. 7, the relative standard deviation   is also considered, assumed to be 0.4, as suggested by 26 for all damage states and building types.Furthermore, there is an uncertainty factor in the response and resistance of the tank (pushover curve).Thus, the relative standard deviation   is also considered, assumed to be 0.25, given the code-designed tank. 26The total uncertainty  tot of the probabilistic seismic demand model 42 is estimated as To verify EFB, a decision model is introduced (B5 in Figure 2).It was decided to use a risk-targeted decision model. 30hus, the derivation begins with the assumption that EFB can be prevented when the probability of the EFB limit state  EFB is less than the acceptable target annual probability of the EFB limit state  EFB, : As the engineers are not as familiar with the probability-based decision model in Equation (17), it is converted to a risk-targeted IM-based decision model, 63 which requires further explanation.First, probabilities can be converted to the mean annual frequency of EFB limit-state exceedance: The decision models in Equations ( 17) and ( 18) are essentially the same because  EFB , when it is low, it is almost equal to  EFB . EFB can be evaluated using the conventional risk equation [27][28][29][30]64 : where H(PGA) is the mean annual frequency of exceedance of PGA; (EFB| = ;  EFB ,  tot ) is the EFB limitstate fragility function, defined by the lognormal cumulative distribution function with two parameters: the median PGA producing the EFB limit state PGA EFB , and the corresponding standard deviation  tot , obtained in the previous step of the proposed methodology.
In the next step in derivation of the risk-targeted IM-based decision model, the target EFB limit-state fragility function (EFB| = ;  EFB, ,  , ) is introduced. 63It is estimated from the known target EFB limit-state risk  EFB, , and the seismic hazard function H(PGA).In this case, Equation ( 19) can be rewritten as The parameters of the target seismic fragility function, PGA EFB,t and  tot, , are not known.However, it can be assumed that  tot, ≈  tot because  tot is not a sensitive parameter with respect to  EFB , and the objective in the design is that  EFB ≈  EFB, .Based on this assumption,  EFB, is the only unknown in Equation (20), and can be calculated iteratively using Equation (20).First, the value of PGA EFB,t must be assumed.Next, Equation ( 20) is evaluated.If the assessed  EFB, is similar to the target EFB limit-state risk, then the assumed value of PGA EFB,t is obtained in the first iteration.Otherwise, a new value of PGA EFB,t is assumed, and calculation of the EFB limit-state probability is repeated.The result of the calculation is the target EFB limit-state fragility function.
Inserting Equations ( 19) and ( 20), both with known parameters, into Equation ( 18), the decision model can be simplified to the risk-targeted IM-based decision model: To fulfil the IM-targeted performance objective from Equation ( 21), the PGA EFB from the estimated EFB limit-state fragility function of the tank must be on the right side of the assessed target EFB limit-state fragility function; in the opposite case, the tank is considered unsafe in the EFB limit state.It is noted that the risk-targeted IM-based decision model was introduced before 63 and has already been successfully implemented for the safety verification of the freeboard of storage tanks. 65However, in this study, the risk-targeted IM-based decision model, Equation (21), is considered for safety verification of the EFB limit state in the wall of storage tanks.

EXAMPLE: UNANCHORED STEEL STORAGE TANK IN ROME
The methodology introduced in Section 2 was applied to an unanchored steel storage tank in Rome, Italy.The site properties (A1) and ground motions at the bedrock level (A2) are presented first.The SRA results (A3) are elaborated.Process B begins by explaining the geometry and material properties of the tank (B1).The loading acting on the tank is defined (B2) and applied to perform pushover analysis (B3).The results of the SRA are coupled with the results of the pushover analysis to obtain the seismic demand model for evaluation of EFB (B4).The EFB limit-state safety is verified by the risk-targeted IM-based decision model (B5).

Site response analysis (Process A)
In this example, it was assumed that the existing unanchored steel storage tank was founded on soil type C, 41 with  = 1.7 t/m 3 and ν = 0.3. 42Knowing the soil density , Equations ( 1) and ( 2) were used to define the reference V s profile of site type C, which was subsequently converted to the interval reference V s profile for 1D SRA (Figure 4).The site interval reference profile was characterised by V s,30 = 325 m/s.The reference shear-wave velocity at bedrock was assumed to be 950 m/s. 42The procedure described in Section 2.1 was used to generate 15 interval V s profiles.Using Equation ( 4), the layer thickness was randomised for 15 profiles by setting the number of layers to 10 and following the boundary criteria for layer thickness described in Section 2.1.The  , profile was randomised using Equation ( 5) for each profile generated from Equation ( 4).The randomly generated  , was converted to V S,Z and to the interval V s profile based on Equation ( 6).The result of this process was a sample of 15 interval V s profiles (Figure 4).As with the soil deposit, randomisation of the bedrock shear-wave velocity  ℎ  was performed using Equation ( 7).The V s,30 of the 15 sample profiles ranged from 320 to 332 m/s (Figure 4); the shear-wave velocity at bedrock ranged from 808 to 1200 m/s.
The ground motions for SRA were extracted from the PEER strong ground-motion database, 36 Engineering Strong Motion (ESM) database 37 and European Strong Motion (EUSM) database. 38The moment magnitude M w interval was constrained between 5 and 8, the epicentral distance between 4 and 65 km, and V s,30 between 650 and 2300 m/s. 42According to these criteria, 47 ground motions were selected from the databases.The PGA of the selected ground motions ranged from 0.01 to 0.53 g (Table 1).With a lack of ground motions recorded on stiff soil and rock, there are little data available for PGA > 0.3 g (Figure 5).Thus, ground motions with PGA between 0.3 and 0.4 g (grey box in Figure 5) were scaled to cover the gap of ground motions with PGA > 0.4 g.Additionally, three samples with PGA between 0.2 and 0.3 g were scaled (green box in Figure 5) to better represent ground motions with PGA between 0.3 and 0.4 g.All considered PGA  values at bedrock are presented in Figure 6.In total, 61 ground motions were considered in the cloud analysis used for SRA.The selection procedure for ground motions for cloud analysis was consistent with the procedure from a previous study. 42he SRA was performed by DEEPSOIL (step A3, Section 2.1) for each equivalent soil profile and ground motion.In total, 15 × 61 site response analyses were performed.The SRA produced the acceleration response spectra S e , corresponding to the tank foundation elevation.In total, 915 acceleration response spectra were obtained at the tank foundation level.The difference between the input acceleration response spectrum at bedrock and the output acceleration response spectrum at the tank foundation level is presented in Figure 7 for the equivalent profile N • 10 and ground motion recorded in Turkey at the Bingol-Bayindirlik Murlugu station (Table 1).
The SRA results were analysed in terms of the relationship between PGA and maximum shear strain throughout the soil profile.In Figure 8, it is observed that by increasing the PGA at the bedrock level, the maximum shear strain of the soil deposit close to the bedrock increases almost linearly with the PGA.However, for PGA greater than approximately 0.25 g, large maximum shear strains were observed at the soil layer close to the surface due to its low initial stiffness G max .Furthermore, the maximum shear strain was always less than 0.35% (Figure 8), which is the maximum shear strain at which the equivalent linear 1D SRA can still be considered acceptable.F I G U R E 9 Quad-linear stress-strain model of S235 steel. 72

Pushover-based seismic performance assessment (Process B)
The pushover-based seismic performance assessment of the tank began with definition of its material and geometric properties (B1, Section 2.2).The radius R and total height H tot of the tank were 27.43 and 15.60 m, respectively.The maximum filling level is defined as 90% of the total height (H = 0.9 H tot ).As the bottom wall thickness t is equal to 0.033 m, the ratios R/t and H/R are equal to 831.27 and 0.51, respectively, indicating that the tank can be classified as a rigid tank according to the working draft of the new Eurocode. 32The tank has a bottom plate thickness t b and annular plate thickness of 0.008 m.Because there is no information about the tank foundation, it was assumed as a concrete ring wall. 66Based on the review of existing literature, 67,68 the concrete ring wall was then considered as the rigid foundation.The tank has wall with varying thickness, as presented in Table 2.The tank material is S235 steel, the most common material for tank construction, as it is economical, and easily fabricated and welded. 52S235 steel is characterised by an elasticity modulus of 210 GPa, a Poisson's ratio of 0.3, a mass density of 7849.7 kg/m 3 and a yield strength of 235 MPa.The quad-linear model of the S235 stress-strain relationship is presented in Figure 9.
Once the tank properties were defined, a 3D non-linear numerical model of the tank was developed in Abaqus (B1, Section 2.2).The mass of the floating roof is only about 0.5% of the total weight.In addition, fluid displacements at the surface are not significantly influenced by the floating roof, 68 because the floating roof is light and flexible. 11,33,69As such, it does not contribute significantly to the increase of the axial stiffness of the wall, which does not apply to the case of the fixed roof. 70As a consequence, the effect of the floating roof on the EFB was neglected.The tank walls and bottom plate were modelled with four-node reduced-integration shell elements (S4R). 16The mesh was characterised by a quad-dominated element shape using an advancing front algorithm. 16The mesh size was 25 cm in the meridional direction of the tank and 70 cm along the tank circumference. 71The stiffening ring at the top of the wall was simulated as 5-mm thickness shell elements S4R.The stiffening ring is of 1114-mm width with tie constraints simulating the truss elements.The uplifting of the tank bottom plate from the assumed rigid foundation was modelled by interface elements.The geometric non-linearity at the tank-foundation interface due to the base uplifting was also considered.The surface-to-surface contact was hard contact, and the friction coefficient between the bottom plate and the rigid foundation was assumed to be 0.3. 16The rigid F I G U R E 1 2 Pushover-based S e -σ relationship, elephant-foot buckling limit-state criterion according to Eurocode and points representing the EFB limit state according to Eurocode and elasto-plastic EFB limit state from pushover analysis.foundation was simulated using an analytical rigid shell.The non-linear material of the tank was defined by the von Mises plasticity model with combined hardening.The tank model is presented in Figure 10.
In step B2, the loading acting on the walls of the tank was defined according to the procedure described in Section 2.2.Based on the configuration of the tank, the patterns of hydrodynamic pressures,  ,ℎ and  , , were evaluated and applied on the tank accounting for the SRSS rule. 32he pushover analysis of the tank (B3) was performed, as described in Section 2.2. Figure 11 shows the von Mises stresses generated by the hydrodynamic pressures acting on the tank at the final step of the pushover analysis.Figure 12 shows the relationship between the IM S e , and the maximum axial compressive stress in the tank wall σ (EDP).The maximum value of σ producing the EFB limit state according to the Eurocode criterion (red dot, Figure 12) occurred at the intersection between the S e -σ curve from the pushover analysis and the Eurocode criterion (Equation 13).The σ producing the EFB limit state was estimated as 36 MPa according to the Eurocode criterion and obtained at S e = 5.40 m/s 2 .However, the σ that caused initiation of plastic strain at EFB region, as observed in the pushover analysis was slightly higher than that F I G U R E 1 3 ln PGA-lnσ points representing cloud-based seismic demand and the linear model obtained from regression analysis in a log-log domain.The intersection of the linear demand model and Eurocode criterion for elephant-foot buckling occurrence defines peak ground acceleration (PGA) for the EFB limit state.estimated according to the Eurocode criterion and was also obtained at a slightly greater spectral acceleration, S e = 6.64 m/s 2 (purple dot, Figure 12).However, the σ producing the EFB limit state according to the Eurocode criterion was considered in this example.Pushover analysis was conducted by increasing the load amplification factor, which is proportional to S e , according to the procedure described in Section 2.2.
It was possible to couple the acceleration response spectra from the cloud-based SRA (A3) and the tank pushover analysis (B3).The SRA provided a seismic demand sample expressed in S e -σ points.S e was evaluated at the ground surface at T ir,h , computed iteratively. 32T ir,h was dependent on the site properties described in Section 3.1, the tank configuration, and the foundation radius, assumed to be 0.5 m greater than the tank radius. 52T ir,h varied slightly for the equivalent profiles; the minimum and maximum T ir,h were 0.129 and 0.134 s, respectively.In total, 15 × 61 S e values were extracted from the SRA (A3, Section 2.1) and used to calculate the corresponding σ values based on the tank pushover analysis (B3, Section 2.2).Note that the 15 × 61 S e values were used to estimate the site amplification factors for the equivalent soil profiles.It was observed that the average amplification factor for low S e values is higher by about 25% compared to short-period site amplification factor F α = 1.6 for soil type C. 41 However, for S e higher than about 0.25 g the site amplification factor decreases.The difference between the maximum values of the site amplification factors is considered acceptable.To some extent, it arises from the Eurocode site amplification factor model that directly accounts for the V s,30 , 41 while the site amplification factors estimated by the SRA account for the V s profiles and other uncertainties related to soil characteristics and ground motions.
In the next step, S e was converted to the corresponding PGA at the rock outcrop (B4, Section 2.2).The outcome of the coupling of SRA and pushover analysis was the sample of the PGA-σ points representing the cloud-based seismic demand (Figure 13).It was then possible to evaluate the logarithm of the median PGA producing EFB limit state  EFB at the intersection of the linear model with the logarithm of the Eurocode criterion for EFB (Figure 13). EFB was computed as 1.14 m/s 2 .The logarithm standard deviation β tot was evaluated according to the procedure described in step B4 (Section 2.2); its value was 0.52, considering the demand variance from the median values   = 0.22, the EFB limit-state definition uncertainty   = 0.40 and the tank capacity uncertainty   = 0.25.
EFB limit-state verification was performed according to step B5 (Section 2.2).The input parameters of the risk-targeted decision model were the target mean annual frequency of EFB limit-state exceedance  EFB, , and the hazard curve at the bedrock at the tank site. EFB, was assumed as 4 × 10 −4 . 73The hazard curve for Rome (Italy) was obtained from the EFEHR open-access database 74 as the arithmetic mean in such a way that considered the effects of the epistemic uncertainties (Figure 14).Based on the mean seismic hazard curve, the adopted  EFB, and the estimated β tot = 0.52, PGA EFB,t F I G U R E 1 4 Seismic hazard curve for bedrock in Rome.

F I G U R E 1 5
Seismic fragility curves for elasto-plastic elephant-foot buckling limit state and Eurocode-based EFB limit state compared to risk-targeted fragility curve.was calculated iteratively using Equation (20) as 5.2 m/s 2 .As PGA EFB was significantly smaller than PGA EFB,t , the EFB limit state was not verified based on Equation (21).Thus, the tank was classified as an under-designed tank with respect to the EFB limit state.Such an outcome is observed graphically in Figure 15, which presents the risk-targeted fragility function and the fragility function of the EFB limit state of the tank according to the Eurocode criterion.The EFB fragility function of the tank appears to the left of the risk-targeted fragility function.
The fragility function of the elasto-plastic EFB limit state of the tank was evaluated (Figure 15).Based on the risk-targeted decision model, either the elasto-plastic EFB limit-state fragility function or the Eurocode criterion EFB limit-state fragility function did not satisfy the safety verification based on Equation (21).However, the two fragility functions are similar, but PGA EFB for the elasto-plastic EFB limit state is higher and amounts to 1.54 m/s 2 .Because the proposed methodology is not limited to using PGA as the seismic IM, the safety verification of the EFB limit state was also performed based on spectral acceleration at bedrock at the tank's impulsive period S e, bedrock .However, changing the ground motion IM for fragility analysis did not change the conclusions resulting from the PGA-based fragility functions.Also in this case, the S e, bedrock,EFB that is equal to 2.6 m/s 2 was significantly smaller than the S e, bedrock, EFB,t = 7.0 m/s 2 .

LIMITATIONS
The proposed methodology focusses on verifying only one failure mode observed in the above-ground liquid storage tanks.Such an approach is common in research, 71 but in the case of practical applications of the seismic performance of the structure, all failure modes (e.g., base uplifting, base sliding, yielding foundation and foundation uplifting) should be verified.
The proposed methodology decouples the problem into the SRA and structural analysis.Consequently, the soilstructure interaction effects are neglected, which is consistent with most of the codes for earthquake-resistant design.However, further research is needed to implement into the proposed methodology the soil-structure interaction effects, especially in the case when the foundation should not be assumed rigid. 75he results of the proposed methodology are as accurate as the pushover analysis can predict the EDPs on the tank.The inertial effects on the liquid and the effects on the evolution of stresses and strains due to time-dependent phenomena, such as fatigue and steel strength due to cyclic loading, cannot be considered in the pushover analysis.However, the present study aims to investigate the elephant-foot buckling phenomenon in the non-linear range.The pushover analysis is thus used in this study because of its computational efficiency compared to the dynamic analysis of the refined 3D non-linear model of the tank.Further studies are needed to better understand how important the hysteretic behaviour of the tank structural elements is with respect to the impulsive response of the tank.
The proposed methodology accounts for the combination of the horizontal and vertical components of the ground motion and the resulting hydrodynamic pressures.In the presented example, the components of the hydrodynamic pressure on the tank due to horizontal and vertical ground motion are assumed to be based on the same level of S e , while total hydrodynamic pressure on the tank is based on the SRSS rule 32 of the two components of the hydrodynamic pressure.However, the proposed methodology allows considering different combinations and orientations of ground motion components that are used to estimate the hydrodynamic pressure for the pushover analysis.
As the study is focussed on EFB, the sloshing phenomenon was neglected.It is assumed that sloshing does not significantly affect the EFB.In the proposed methodology, only the impulsive component of tank response contributes to the EFB.The variation of filling level over time that affects the impulsive response component of the liquid was also disregarded.For simplicity, the maximum filling level was considered in this study, which refers to 90% of the tank wall height. 7,53

CONCLUSIONS
A methodology was presented for risk-targeted seismic performance assessment of EFB of liquid storage tanks.It consists of cloud-based 1D site response analyses, pushover analysis using a refined 3D non-linear model and the risk-targeted IM-based decision model.The SRA couples the seismic hazard analysis for the bedrock with the spectral acceleration at the impulsive period of the tank and the tank foundation elevation.The spectral acceleration is then used to define the hydrodynamic pressure on the tank as defined in the working draft of the new Eurocode.Finally, the seismic demand on the tank wall is calculated by the pushover analysis using the refined non-linear structural model of the tank.The proposed methodology can be used for risk assessment of EFB in liquid storage tanks, but this study involves a more simplistic risk-targeted IM-based decision model.To verify safety from EFB , the engineer must verify that the seismic IM causing EFB is higher than the seismic IM corresponding to the target seismic risk.
In the proposed methodology, the axial stresses of the tank wall are obtained from the pushover analysis that is related to the spectral acceleration at the tank's impulsive period.The elasto-plastic EFB is simulated in the pushover analysis, while the conversion between the spectral acceleration at the tank's impulsive period and the seismic IM used in the seismic hazard analysis is obtained from the cloud-based SRA.
The investigation of other failure modes of steel storage tanks due to seismic action (e.g., base uplifting, base sliding and top wall buckling) was beyond the scope of this study.Additional research is needed to investigate if the proposed methodology can be extended to the verification of other failure modes of the tanks.
An example demonstrated that sampling of the interval shear-wave velocity profiles and 1D SRA are relatively simple to implement if SRA computations are automated.For the investigated tank it was demonstrated that the EFB limit state defined in the Eurocode is only slightly conservative concerning the occurrence of elasto-plastic EFB.

A C K N O W L E D G E M E N T S
The work presented herein was carried out with a financial grant from the H2020 Marie Skłodowska Curie Innovative Training Network as part of the research project, XP-RESILIENCE (G.A. 721816), and the Slovenian Research Agency, as part of the Earthquake Engineering Research Program (P2-0185).This support is gratefully acknowledged.The authors also thank Assoc.Prof. Primož Može for discussions on the modelling and analysis performed in Abaqus.

F I G U R E 4
Interval reference V s profile and equivalent generated profiles.

F I G U R E 5
Sample of peak ground acceleration (PGA) values of 47 ground motions selected from strong ground-motion databases (European Strong Motion [EUSM], Engineering Strong Motion [ESM] and PEER).

F I G U R E 6
Sample of peak ground acceleration (PGA) values of 61 ground motions selected from three databases and scaled ground motions used in 1D site response analysis (SRA).F I G U R E 7 Input acceleration response spectrum at bedrock and output acceleration response spectrum at tank foundation level for equivalent profile N • 10 and ground motion recorded in Turkey at Bingol-Bayindirlik Murlugu station.F I G U R E 8 Maximum shear strain as a function of depth and peak ground acceleration (PGA) for all equivalent soil profiles and ground motions in cloud-based site response analysis (SRA).TA B L E 2 Thickness of tank wall as a function of elevation.

F I G U R E 1 0
3D non-linear numerical model of the investigated tank.F I G U R E 1 1Von Mises stresses generated by the hydrodynamic pressures acting on the tank at the final step of the pushover analysis.
Selected ground motions for SRA and their basic characteristics.
TA B L E 1PGA, peak ground acceleration; SRA, site response analysis.