Energy‐compatible modulating functions for the stochastic generation of fully non‐stationary artificial accelerograms and their effects on seismic site response analysis

The reliability of non‐linear dynamic analysis aimed to predict the seismic performance of structural and geotechnical systems, as well as their dynamic interaction, is significantly affected by the selection of suitable input motions. Several procedures for selecting actual earthquake records compatible with the seismic hazard at the site of interest and for evaluating synthetic accelerograms consistent with a seismological framework are available in the literature. However, the degree of uncertainty affecting these approaches might lead to unreliable performance predictions especially for strongly non‐linear problems, such as those involving the response of soils to cyclic and dynamic loading. In this vein, the paper presents a procedure for generating fully non‐stationary ground motions ensuring energy compatibility with a target motion. Two different approaches have been introduced: i) the time window method based on the central finite difference approach and ii) the simplified intensity‐compatible approach in which a closed form expression is provided to evaluate a modulating function that depends on the Arias intensity and on the strong motion duration of a target motion. To highlight the reliability and the accuracy of the proposed procedure, several sets of spectrum‐compatible artificial accelerograms, generated starting from a rock outcropping motion, were used as input motions in a series of one‐dimensional site response analyses. The analysis results are presented and discussed in the paper highlighting the influence of the main features of the proposed generation procedure on the variability of the computed site response.


NOVELTY
• A model of the evolutionary power spectral density function that depends only on the frequency of peaks, zero level up crossings and on the total energy of a target accelerogram, has been used to generate artificial accelerograms consistent with actual seismic record. • A novel energy compatible modulating function has been calibrated by using: i) the time window method (based on the central finite difference approach); ii) the simplified intensity-compatible approach that provides a closed form expression depending on the Arias intensity and on the strong motion duration of a target motion. • Artificial accelerograms are used as input motion in 1D site response analyses.

INTRODUCTION
In most of seismic codes, the earthquake-induced ground shaking is generally represented in terms of pseudo-acceleration response spectrum. However, for strongly non-linear problems, such as those involving the response of soils to cyclic and dynamic loading, the use of the response spectrum technique may be inadequate 1,2 and non-linear dynamic analyses are required, implying the need of selecting suitable input motions. The definition of proper sets of acceleration time histories is a crucial step since it significantly affects the reliability of the analyses aimed to predict the seismic performance of structural and geotechnical systems, as well as their dynamic interaction. In this framework, the increasing availability of strong motion records makes the use of actual accelerograms, recorded during large earthquakes, an attractive option to define the input motions. The use of properly selected sets of acceleration time histories allows accounting for the uncertainties involved in the assessment of the seismic hazard at the site of interest. Therefore, to account for the inherent variability of the expected ground motion, large sets of recorded accelerograms should be selected and used in the seismic analyses carried out with predictive purposes. [3][4][5] The assessment of the safety conditions of seismically excited structures has been addressed by Muscolino et al. 3 considering the inherent uncertainties influencing the definition of the ground motion accelerations. Specifically, two large sets of accelerograms, recorded con rock subsoil, have been analysed by these Authors 3 to predict the upper and lower bounds of the reliability function of structural systems subjected to earthquake accelerograms.
The results of the procedures devoted to the selection of actual accelerograms are influenced by multiple sources of uncertainties mainly related to the definition of the seismic hazard at the site of interest and to the effect of non-linear soil behaviour influencing the site response and the expected motion at the ground surface or at the rock outcrop. Furthermore, the introduction of proper compatibility criteria of the selected records with the frequency and energy content of an expected motion (assumed as target in the selection procedure) is also required. Since the target ground motion is usually defined in terms of code-prescribed elastic response spectra, the geotechnical properties of soil affecting the nonlinear site response should be accounted for in the definition of the interval of vibration periods adopted for checking the compatibility of the selected records. 6,7 Different selection procedures have been proposed in the literature either with structural or geotechnical purposes. [8][9][10][11][12][13][14][15][16][17][18] However, depending on the adopted compatibility criteria and on the characteristics of the target motion, the selection of a suitable number of compatible actual accelerograms might result impossible without applying large scale factors which distort the characteristics of the selected motions. 14,19,20 This may lead to unrealistic inputs and hence to unreliable predictions of the seismic performance of structural and geotechnical systems 14,[19][20][21][22][23] ; for the latter case, Biondi et al. 21 and Aliberti et al. 23 have shown that the influence of the input motion selection criteria on the prediction reliability is relevant regardless the complexity level of the analyses.
Due to the nature of earthquakes, the rational way to assess the performance of earthquake-resistant structures consists in the use of the stochastic approach for modelling the dynamic action. The stochastic ground-motion models frequently used by engineers and seismologists include both physics-based [24][25][26][27] and record-based 28-31 models. Therefore, alternatively to actual records, sets of synthetic accelerograms obtained through seismologically based numerical models and accounting for site-to-source effects could be used as input motions. However, a proper selection of the source parameters (generally based on data related to previous earthquakes) invariably conveys a high degree of uncertainty; also, the ground motions predicted for future earthquake scenarios can be highly sensitive to the specification of these parameters and require an expert engineering judgment. 17,18 For these reasons, the use of synthetic accelerograms is infrequent in the dynamic analyses carried out with engineering purposes. Artificial accelerograms represent a valuable alternative to the previously described approaches and can be suitably adopted for several geotechnical and structural issues involving large non-linear behaviour of soils and structural materials. Traditionally, artificial accelerograms are generated using numerical procedures aimed to match a target response spectrum for which a power spectral density (PSD) function is preliminarily evaluated and accelerogram samples are then derived as sinusoidal signals having random phase angles 32,33 ; alternatively, the S-transform 34 and the wavelet analysis [35][36][37][38][39][40][41][42] can be also adopted to generate artificial accelerograms. Genovese et al. 42 used the circular wavelet transform to decompose a recorded accelerogram into the superposition of complex-valued harmonic wavelets with complex-valued combination coefficients, which are then randomized through a generalization of the well-known Shinozuka's formula. 43 To generate artificial accelerograms, Muscolino et al. 44 proposed a new model of the evolutionary power spectral density (EPSD) function of zero mean Gaussian processes that depends only on the frequency of peaks, zero level up crossings and on the total energy of a target accelerogram.
Artificial accelerograms generated by applying most of the procedures available in the literature are not able to capture the evolution of the time axis up-crossing which clearly reflects the evolutionary frequency content of the ground motion recorded in sites having different geotechnical properties. Moreover, it has been recognized that artificial accelerograms generated by applying stationary models (extensively suggested by several international seismic codes) are characterized by an excessive number of cycles and consequently possess unreasonably much higher energy content with respect to actual ones. 18 Accordingly, these kind of artificial accelerograms should not be used in the dynamic analyses aimed to predict the response of soil deposits and geotechnical systems which are generally governed by large non-linearity effects and frequently involve large plastic strains.
In this vein, a novel procedure for generating samples of fully non-stationary zero-mean Gaussian processes is described in this paper. Starting from a properly selected target accelerogram, the proposed procedure allows detecting a fully nonstationary model of the ground motion for which the target accelerogram may be considered as one of its samples. In the proposed procedure, the EPSD-based model proposed by Muscolino et al. 44 has been improved, by a proper calibration of a novel energy compatible modulating function, and two different methods have been introduced: i) the time window method (TWM) in which the central finite difference approach has been used to capture the temporal evolution of the absolute amplitude of the selected accelerogram; ii) the simplified intensity-compatible approach (SIA) in which a closed form expression is given for the evaluation of a modulating function that depends on the Arias intensity and on the strong motion duration of the target motion. Unlike available studies, the time and frequency characteristics of the proposed ground shaking model are calibrated separately. This peculiarity allows to easily calibrate the modulating function in the time domain by using the proposed TWM or SIA. Furthermore, the closed form expression given by SIA allows to evaluate a modulating function that depends only on seismic parameters (easily predictable for a given site by using empirical attenuation relationships 45 ) instead of using an acceleration time history. The proposed modulating functions make easier the application of the EPSD model presented by Muscolino et al. 44 given that the proposed expressions (TWM and SIA) avoid solving an optimization problem for the calculation of the various parameters on which depends the modulating function of the aforementioned model. 44 Furthermore, differently from the procedure proposed by Muscolino et al., 44 an additional spectrum-compatibilization procedure for the EPSD-based approach has been introduced in the present study. Specifically, an iterative procedure is applied to obtain the compatibility between the mean elastic response spectrum of the set of generated accelerograms and the spectrum of the target accelerogram.
The key aspects of the model proposed by Muscolino et al. 44 have been reported in Sections 2 and 2.1 to make the present paper self-contained. The proposed time window method (TWM) and simplified intensity-compatible approach (SIA) are described in Section 2.2, while the procedure adopted to generate artificial accelerograms is described in Section 2.3. To the authors knowledge, the use of sets of fully non-stationary generated accelerograms in the dynamic analyses devoted to the prediction the seismic response of geotechnical systems is not frequent. Accordingly, the reliability and accuracy of the proposed procedure is illustrated in the Section 3 of paper with reference to one-dimensional seismic site response analyses carried out accounting for non-linear soil behaviour. The results are presented and discussed to highlight the influence of the proposed generation procedure on the computed profiles of maximum acceleration, stratigraphic amplification factors and earthquake-induced shear strains and on the main characteristics of the motion predicted at the ground surface. The analysis results show that the proposed procedure is able to generate artificial accelerograms that satisfactorily reproduce both the time-varying amplitude and frequency content of the target accelerogram, leading to a reliable prediction of the site seismic response.

FULLY NON-STATIONARY ENERGY-COMPATIBLE ARTIFICIAL ACCELEROGRAMS
Earthquake-induced ground accelerations can be modelled as samples of a zero-mean fully non-stationary Gaussian random process. This can be evaluated as the product of a sample of a zero-mean stationary Gaussian process by a suitably selected time-frequency dependent modulating function. A useful approach to generate samples of fully non-stationary Gaussian processes is based on the evolutionary spectral representation, that requires the introduction of the Evolutionary Power Spectral Density (EPSD) function. [46][47][48] In this paper, the formulation recently proposed by Muscolino et al. 44 is adopted to this purpose. According to this formulation, an actual (recorded) accelerogram̈( ), of overall duration D, can be considered as one of the samples of a fully non-stationary process 0 ( ) defined as the sum of zero-mean Gaussian uniformly modulated random processes F 0,k (t). Each of these processes consists of the product of a non-negative deterministic modulating function a(t) by a stationary zero-mean Gaussian filtered sub-process, X k (t): In Equation (1), the process 0 ( ) is obtained by dividing its overall duration D into n contiguous time intervals of amplitude ΔT k = t k -t k-1 and W(t) is a window function, which is unitary for t k-1 ≤ t< t k and null elsewhere.
The one-sided EPSD function of the zero-mean Gaussian fully non-stationary process 0 ( ) is given by: where k ( , ) and k ( ) are the one-sided EPSD and PSD functions of the sub-process X k (t), respectively.
In the time interval [t k-1 , t k ), the k th sub-process X k (t) is characterized by the following one-sided PSD function: where ω L,k and ω H,k represent the k th frequency control of the second order low-pass and first order high-pass Butterworth filters, respectively, and The parameters ρ k and Ω k introduced in Equation (4) represent measures of the frequency bandwidth and of the predominant circular frequency of the k th filtered stationary process, respectively. The coefficient β k introduced in Equation (3) ensures that the sub-process X k (t) possesses a unitary variance, that is 2 k = 1, and can be evaluated as: The coefficients appearing in the closed form expression of Equation (5) can be calculated as: The so defined fully non-stationary stochastic process F 0 (t) is able to capture simultaneously the time-varying intensity and the time-varying frequency content of the target accelerogram̈( ), although its single processes are individually uniformly modulated. The similarity between̈( ), and each sample of the process F 0 (t) can be improved by separately fine-tuning the modulating function and the frequency content of F 0 (t). Specifically, since each of the subprocess X k (t) introduced in Equation (2) is characterized by a one-sided PSD function, k ( ) with unit variance, the modulating function and the main parameters characterizing each PSD function can be estimated separately; this feature represents the main novelty of this model in comparison with existing studies.

Parameters of the sub-process PSD function
Once the overall duration D of the target accelerogram is divided in n contiguous time intervals ΔT k , the one-sided PSD function X k ( ) of each stationary sub-process X k (t) can be achieved by capturing, in each of these intervals, a group of waves possessing the specific frequency distribution of the target accelerogram in the same time interval. In this vein, the spectral parameters ρ k , Ω k , ω H,k and ω L,k introduced in Equations (3)-(5), must be properly evaluated. According to Muscolino et al., 44 the control frequencies of the k th pair of Butterworth filters in the time interval [ −1 , ) are given by H,k = 0.1 Ω k and L,k = Ω k +0.8 k , and the predominant circular frequency and the circular frequency bandwidth can be evaluated as: where + 0,k and k are the number of zero-level up-crossings (i.e. crossings of the time axis with positive slope) and the number of peaks of the target accelerogram in the time interval [t k-1 , t k ).

Modulating function: energy-compatible approach
To estimate the modulating function a(t), the current value of the cumulative intensity I 0 (t), proportional to the cumulative energy associated to the target accelerogram, is firstly evaluated: together with its cumulative value I 0 obtained setting t = D in Equation (9).
Owing to the unit variance of the stationary sub-process X k (t), the expected cumulative energy of the fully non-stationary stochastic process F 0 (t) can be evaluated as: In Equation (10), the linear operator E⟨⋅⟩ provides the expected value of the quantity within angle brackets. The expected cumulative intensity of the stochastic process F 0 (t) can then be calibrated against the deterministic function I 0 (t). To this purpose a least-square fitting can be adopted even if it can be relatively cumbersome for practical engineering applications. 44 Thus, two alternative methods will be introduced in the following Sections.

2.2.1
Time window method (TWM) The first method, here referred as time window method (TWM), evaluates the modulating function analysing the time evolution of the absolute amplitude of selected recorded accelerogram, using the following relationship: where I m (t) is the cumulative intensity of the modulating function a(t). The relationship in Equation (11) holds because all the stationary sub-processes X k (t) have unit variance. Thus, the modulating function a(t) can be evaluated analysing the absolute amplitude of the target accelerogram and the following relationship holds: Since I 0 (t) is a monotonic increasing function, using the central finite difference approach, Equation (12) reduces to: where t j = j Δt (with j = 0, . . . , N) and the time window Δ̃> Δ must be opportunely selected, Δ̃being a multiple of the sampling step Δt of the target accelerogram while N is the number of discrete time instants of the target signal. To avoid inconsistent results, must also be satisfied the following conditions a(t 0 ) = a(D) = 0 and Δ̃≤ j = j Δ ≤ − Δ̃.

Simplified intensity-compatible approach (SIA)
To avoid the numerically evaluation of time derivatives involved in the previously described TWM, a further method, herein denoted as simplified intensity-compatible approach (SIA), is proposed requiring the a priori definition of a modulating function. Herein, the following Jennings-type modulating function 50 is adopted: with: The dimensionless coefficient 0 , introduced in Equation (14), can be evaluated by substituting this equation into Equation (11) and equating the result to the expected cumulative intensity I 0 of the target accelerogram (obtained setting t = D in Equation (11)). Furthermore, since, the integral of the function J 2 (t, t 1 , t 2 , D) can be evaluated in the following closed form: the dimensionless coefficient 0 , taking into account that exp (−20) ≪ 1, can be evaluated as: Introducing the Arias intensity 51 of the target accelerogram I A = I 0 π/(2 g), and assuming t 1 = t 05 and t 2 = t 95 Equation (17) leads to: where g is the gravity acceleration, and t 05 and t 95 being the time instants defining the strong motion duration 52,53 D 5-95 of the target accelerogram, for which the 5% and the 95% of the total intensity I 0 has been released, respectively. According to this formulation, through I A , t 05 , t 95 , D, and D 5-95 , the main characteristics of the target accelerogram (hereafter denoted as target motion) could be accounted for in generating the samples (i.e. the so-called artificial accelerograms) of the fullynon stationary zero-mean Gaussian process. It is worth noting that for these quantities several empirical attenuation relationships have been proposed in the literature allowing to reliably predict the values of these parameters for the ground motion expected at a given site. 45

Fully non-stationary samples
Once the parameters characterizing the fully non-stationary zero-mean Gaussian process F 0 (t) have been estimated, its samples can be generated in such a way that the target motion̈( ) can be considered as one of them. By using the Shinozuka's formula, 43,54,55 the i th sample of the stochastic process F 0 (t), can be evaluated as: assuming a frequency step Δω = 2π/D, an upper cut-off circular frequency ω N = π/Δt and the random phase angles, ( ) , uniformly distributed over the interval [0 − 2 ). It follows that m N = ω N /Δω and Δt = π/(4ω N ). 56,57 In current practice the seismic actions to be used in seismic design are generally provided through the displacement ( ) d or pseudo-acceleration ( ) pa elastic response spectra. The generation of a set of fully non-stationary accelerograms compatible with the response spectrum of the target motion can be evaluated according to the procedure described in Genovese et al. 42

SITE RESPONSE ANALYSIS USING ARTIFICIAL INPUT MOTIONS
To highlight the reliability and the accuracy of the proposed procedure in generating sets of artificial accelerograms with characteristics similar to those of the target one, envisaging responses with a reduced variability, an application to onedimensional (1D) equivalent-linear site response analysis is presented and discussed herein. Starting from a target motion, the proposed procedure was used to generate different set of artificial accelerograms which have been used as input motions in a set of seismic site response analyses carried out with reference to an ideal soil profile. In the following Sections, the main features of the selected soil profile and target motion are preliminary introduced. Then, the two proposed energy compatible approaches were applied generating sets of accelerograms using both the time window method (TWM) and the simplified intensity compatible approach (SIA). Once the iterative procedure to obtain the spectrum compatibility was applied, the accelerograms of the corresponding final sets have been used as outcropping input motions in the site response analyses. The analysis results are presented and discussed in terms of variability of predicted site responses in comparison with those computed using the target motion.

Soil profile and numerical modelling
To focus only on the characteristics of the artificially generated input motions and to point out their influence on the predicted site response, only one soil profile was considered in the analysis. It is a simple ideal soil profile schematizing the actual conditions of an actual site characterized, in the shallowest portion, by coastal deposits consisting of sand and gravel with little or no fine content. Starting from the results of a set of down-hole and MASW measurements carried out in the area of interest, an ideal profile of the shear wave velocity was defined fitting the whole set of experimental data through the following expression 58 : where V S0 is the value of the shear wave velocity V S at the ground surface (z = 0), H = 30 m is the thickness of the soil column considered in the site response analysis and the parameters m and α allow describing the degree of heterogeneity of the considered V S profile. The fitting procedure leads to the values V S0 = 150 m/s, m = 0.75 and α = 7.59 which allow averaging the actual sets of V s profiles. The considered soil profile is characterized by an average shear wave velocity in the topmost 30 m equal to about V s,30 = 400 m/s and, thus, belongs to the soil class B (360 m/s < V s,30 < 800 m/s), according to EC8. 59 The first three natural periods of horizontal elastic vibration, evaluated starting from the largest peaks of the undamped elastic transfer function computed between the top and the base of the considered soil column, are equal to about T 1 = 0.29 s, T 2 = 0.13 s and T 3 = 0.08 s, respectively, corresponding to the vibration frequencies f 1 = 3.47 Hz, f 2 = 7.78 Hz, f 3 = 12.57 Hz. The site response analyses were carried out using the computer code DEEPSOIL vs. 6.1. 60 Since the analyses aims only to examine the influence of the characteristics of the generated artificial accelerograms on the computed site response, the simplified linear-equivalent approach was used, and some simplified assumption listed below were introduced in the definition of the geotechnical model. Due to the minimum value V S0 = 150 m/s of the shear wave velocity in the considered soil profile and to the largest frequency of the selected target motion (equal to about 40 Hz), the 1D numerical model was discretized into 0.5 m thick layers to avoid numerical distortion of the propagating waves. This discretization ensured that the maximum height of each element of the numerical model was smaller than 1/6 of the wavelength associated to the highest frequency of the considered input motions. 61 Starting from the results of the site investigations, an average value of the soil unit weight γ = 19 kN/m 3 together with a dry soil condition were considered.
To avoid the influence of other parameters affecting the site response and, again, to focus only on the influence of the characteristics of the generated artificial accelerograms on the computed site response, a rigid bedrock was assumed. Also, the effective shear strain ratio was evaluated starting from the moment magnitude of seismic event considered in the selection of the target motion and, according to Kramer, 57 a frequency-independent complex shear modulus was used. Finally, the curves proposed by Vucetic and Dobry 62 for a plasticity index PI = 0 have been adopted to describe the shear modulus and damping ratio variation with the strain level.

Target motion
The acceleration record̈( ), assumed as target motion is shown in Figure 1A.  Figure 1B and 1C, respectively. The Fourier amplitude spectrum (grey line) of the target motion together with the smoothed one (black line) are shown in Figure 1D. The vertical lines in Figure 1D represent the first three vibration frequencies f 1 , f 2, f 3 of the soil deposit. The predominant period, 57 corresponding to the largest peak of the smoothed spectrum, is T p = 0.615 s (f p = 1/T p = 1.62 Hz). The mean period 64 is equal to about T m = 0.37 s.

Generation of input motions
According to Section 3, the application of the EPSD method requires that each subprocess ( ) is characterized by a one-sided PSD function k ( ) with unit variance. Accordingly, the modulating function ( ) and the main parameters characterizing the PSD function have been estimated separately. F I G U R E 2 Target motion: one-sided PSD function.

Evaluation of the PSD function parameters
According to Muscolino et al., 44 the PSD function ( ) associated to the Gaussian sub-process, X k (t) depends only on the number of peaks P k and zero-level up crossings + 0,k that occur in the k th time interval ΔT k in which the target accelerogram has been subdivided. Muscolino et al. 3 showed that most of the accelerograms recorded on rock outcropping sites are characterized by an almost linear variation of the cumulative number of zero-level up-crossings within the strongmotion duration D  . Furthermore, [65][66][67] pointed out that the zero-level up crossings + 0 and the frequency content of a given accelerogram are strongly correlated. Since accelerograms recorded on rock sites are also characterized by an almost constant frequency content, they can be modelled as samples of zero-mean uniformly modulated or quasi-stationary Gaussian processes. 3 As a consequence, the selected motion can be characterized by only one PSD function ( ) = 1 ( ) (i.e., a single sub-process X 1 (t) was used, meaning that the random process F 0 (t) is uniformly modulated). Thus, the process 0 ( ), of duration D, was obtained considering only one time interval equal to the overall duration (ΔT k = D). The one-sided PSD function 1 ( ) computed for the target motion is shown in Figure 2, while Table 1 lists the parameters useful for its description (Equation (3)).
Notice that this procedure can also be applied to accelerograms recorded on soft soils and to near-fault pulse-like motions. It should also be noted that for soft soils the recorded accelerograms are often characterized by a non-linear TA B L E 1 Main parameters of the one-sided PSD function computed for n = 1 and ΔT 1 =D. trend of the cumulative zero-level up-crossing function and by a non-constant frequency content. In this case, to obtain accurate results, it is necessary to subdivide the target accelerogram into several different time intervals (each containing a number of zero level up crossings equal to at least one) and then, according to Equations (2)-(8), to evaluate a PSD function for each of the time interval into which the target motion has been subdivided (see, e.g., Genovese et al. 68 ). The same considerations also apply for near-fault pulse-like motions recorded on soft soils.

Energy-compatible modulating functions
In the Energy compatible methods described in Section 2, the temporal variation of the amplitudes of the generated samples is obtained through a proper estimation of the modulating function a(t) which in turn depends on the cumulative intensity 0 ( ) of the target motion according to Equations (11)- (13) and Equations (14)- (18). Both the proposed TWM and SIA were adopted herein to this purpose. Using the TWM, eight values of the time window Δ̃, ranging from 0.25 s to 4.0 s, have been used and the corresponding modulating functions are shown in Figure 3 (solid black lines) together with the absolute value of the target motion (grey line); in the same figures, the plots of the integral of the square of the modulating function (herein referred as cumulative intensity of modulating function I m (t)) are superimposed to the trend of the cumulative intensity I 0 (t) of the target motion. It can be observed that the cumulative functions evaluated by using the smallest time windows (Δ̃≤ 1.5 s) are generally characterized by trends very close to the cumulative intensity of the target motion; some differences in the trends can be detected for Δ̃larger than 2 s. However, in all the considered cases, the final values of the cumulative intensity functions coincide with the total intensity of the target motion, implying the same areas under the modulating function plots and, thus, the same total energy. The modulating function evaluated through the SIA is shown in Figure 4. Since for the target motion it is t 1 = t 05 = 5.16 s and t 2 = t 95 = 12.48 s, Equation (18) leads to α 0 = 0.45 m/s 2 . In this case, despite the final values of I 0 (t) and I m (t) at t = D are practically coincident, the differences in the trends are more relevant in comparison with those computed using the TWM.
The EPSD functions (see Equation (2)) computed using the TWM and the SIA are plotted in Figures 5 and 6, respectively.

Fully non-stationary samples
Once the parameters characterizing the fully non-stationary zero-mean Gaussian process, F 0 (t), are estimated, it is easy to generate its samples in such a way that the selected accelerogram can be considered as one of its own samples. By using Equation (19), a set of N s = 100 samples has been generated for each of the nine modulating functions previously introduced (eight with reference to the TWM and one for the SIA) by using a frequency step Δω = 0.157 rad/s, a cut-off circular frequency (equal to the Nyquist's frequency 57 ) ω N = 157 rad/s, and a number of harmonic terms in the right-hand side of Equation (19) Spurious baseline trends, usually noticeable in the displacement time-history obtained by double time-integration of acceleration records, have been removed by applying the baseline correction. Specifically, a best-fit polynomial curve of order 3 is determined for each child sample through a least-squares regression analysis and is then subtracted from the acceleration time history. 42 Figures 7 and 8 show the comparison between the cumulative intensity function of the target accelerogram I 0 (t) (red line) and those of the accelerograms generated using the TWM (Figure 7) and the SIA (Figure 8), respectively. In both these figures, the grey area defines the envelope of the cumulative intensity functions computed for each set of N s artificial accelerograms; the corresponding mean trends with their confidence intervals (computed as the mean plus/minus one standard deviation σ of the data) are also represented with a continuous and a dashed black line, respectively. In the time axes of all the plots, the time interval corresponding to the strong motion duration D 5-95 of the target motion is also identified by the grey vertical lines.

F I G U R E 3
Comparison between the absolute values (grey line) of the target motion, its cumulative intensity I 0 (t), the modulating function a (t) and its cumulative intensity I m (t) computed using the TWM with different time windows. From Figure 7 (TWM), it is apparent that, within the strong motion duration, a general very good agreement is obtained between the cumulative intensity of the target accelerogram (hereafter referred as target intensity) and the mean trend of the intensity functions computed for the N s artificial accelerograms; some differences can be observed only for the larger values of the time window (Δ̃= 3 and 4 s). Regardless the values of the time window, the final (cumulative) values of the computed mean trends are quite coincident with the cumulative target intensity I 0 .
Concerning the SIA (Figure 8), despite the final values of the intensity of the target and mean functions are again coincident, more relevant differences in the trends can be observed within the strong motion duration. Finally, a comparison of the data plotted in Figures 7 and 8 allows to verify that the confidence intervals (and thus the standard deviation) of each set of accelerograms are slightly influenced by the adopted procedure (TWM or SIA) and that the trend of the mean F I G U R E 4 Comparison between the absolute values of the target accelerogram (grey line), its cumulative intensity I 0 (t), the modulating function a (t) and its cumulative intensity I m (t) computed using the SIA.
cumulative intensity function evaluated using the SIA is close to that obtained using the TWM with the largest time window (Δ̃= 4 s).
The comparison between the target function and the mean trends computed for set of N s artificial accelerograms generated using the TWM and the SIA is also presented in Figures 9 (TWM) and 10 (SIA) in terms of cumulative zero-level up-crossings functions. Regardless the considered approach, in both these figures it is apparent an almost linear trend of the plots describing the mean valuē+ 0 ( ) and the corresponding confidence intervals (̄+ 0 ( ) ± ) of the zero-level upcrossings computed for each set of artificial accelerograms. These linear trends are consistent with the assumption of a single time interval in Equation (19), meaning that the same expected zero-level up-crossing rate is assumed for the whole duration of the generated accelerograms. By comparing Figures 9 and 10, it can be observed that, within the strong motion duration, regardless the adopted approach (TWM or SIA), there is a very good agreement between the trends of the zerolevel up-crossings computed for the target motion and for each set of the artificial accelerograms. It is worth nothing that the differences between the trend of the mean valuē+ 0 ( ) and the trend corresponding to the target motion should be ascribed only to the counting of zero crossings in the time intervals preceding (t < t 05 ) and following (t > t 95 ) the strong motion duration.

Spectrum-compatible artificial accelerograms
The procedure used in Section 3.3.3 returns samples of a fully non-stationary random process, such that, in statistical sense, their cumulative intensity functions and zero-level up-crossing functions match those of the target accelerogram. However, the aforementioned functions are not always sufficient to satisfactorily characterize the dynamic action for engineering applications. In current practice, the response spectrum is widely utilized to characterize the seismic action for design purposes. For an actual record the pseudo-acceleration response spectrum ( ) pa is computed, at various natural frequencies ω r or periods T r = 2π/ω r , multiplying 2 r by the absolute maxima of the dynamic response computed for single degree of freedom oscillators having unitary mass and assigned damping ratio 0 .
Conversely, the pseudo-acceleration response spectrum of artificial accelerograms S ( 0 ) pa , with EPSD function 0 defined by Equation (1), is evaluated as the arithmetic mean of absolute values of the responses multiplied by 2 r : In Equation (21), N s is the number of artificial accelerograms generated by Equation (19) and According to the proposed procedure, the set of generated accelerograms is characterized, on average, by an energy and a frequency content very close to those of the target motion. Furthermore, on average, the cumulative number of the zero up-crossings of the time axis (strictly connected to the frequency evolution of the process) and the cumulative intensity (strictly connected with the time evolution of the process) of the target motion are also suitably reproduced. Despite these characteristics, for each of the generated sample it is: meaning that the mean spectra of the generated accelerograms do not satisfy the compatibility condition with the displacement and pseudo-acceleration response spectra of the target motion (hereafter referred as target spectra).
The gap between the target spectrum ( ) pa ( r , 0 ) and the mean elastic pseudo-acceleration response spectrum of the set of artificial accelerograms pa ( r , 0 ) are the value, at the (j-1) th iteration, of the PSD function and the elastic pseudo-acceleration response spectrum of each sub-process X k (t), respectively.
In Equation (25), the following conditions are assumed at the first iteration (j = 1): At the j th iteration, the i th sample of the non-stationary process can be evaluated as: and, afterwards, baseline adjusted. [70][71][72] Obviously, at the j th iteration, the fully non-stationary process is characterized by the following EPSD function: F I G U R E 7 Cumulative intensity functions of the target motion and of the initial set of artificial accelerograms generated using the TWM, with different time windows together with the indication of the target D 5-95 delimited by two vertical lines. and cannot match, at the same time, the energy and the frequency content of target motion as well as the spectrum compatibility to its elastic response spectrum. However, the proposed procedure provides a good compromise in satisfying both these three compatibility issues. A set of preliminary numerical analyses (not presented herein due to lack of space) showed that for the proposed iterative procedure few iterations are required to attain, at some frequency (or period) intervals and within an acceptable tolerance, an equivalence between the target spectrum and the mean elastic response spectrum of the set of generated artificial accelerograms by Equation (28). In particular, to satisfy the spectrum compatibility between the mean elastic response spectrum of each set of artificial accelerograms, pa ( r , 0 ), and the elastic response spectrum of the target motion, ( ) pa ( r , 0 ), the procedure described above was applied with only two iterations. Accordingly, different sets of spectrum-adjusted accelerograms have been generated through Equation (28); specifically, Specifically, with reference to the initial sets of artificial accelerograms (hereafter referred as 0 th iteration), Figure 11A shows a comparison between the elastic response spectrum of the target motion (hereafter referred as target spectrum) and the average spectra of the initial sets of the N s artificial accelerograms; the corresponding relative difference is plotted in Figure 11B versus the vibration period. As it can be observed, regardless the value of the time window considered in the generation procedure, relevant differences between the target and the mean spectra are apparent with relative differences often larger than about ± 30% regardless the vibration period. The results obtained in terms of spectrum compatibility after the two iterations involved in the proposed procedure are presented in the Figures 11 C,D (first iteration) and 11 E,F (second iteration). The influence of the considered time window becomes almost negligible after the first iteration and, after the second iteration, a very good agreement is obtained between the average spectra computed for each set of artificial accelerograms and the target spectrum; specifically, with few exceptions, the relative difference is always smaller than about ± 15% regardless the vibration period. The influence of the time window, which appear negligible with reference to the spectrum compatibility check, is relevant if the set of generated accelerograms and the target motion are compared in terms of amplitude, energy and frequency content.
The results obtained using the SIA are presented in Figure 12A where the target spectrum is compared with the average spectra of the initial set (0 th iteration) of the artificial accelerograms and with the two sets obtained after the first and the second iterations; the corresponding percentage relative difference is plotted in Figure 12B.
A very good agreement between the target and the computed mean spectra is apparent for both the first and the second iterations with values of the percentage relative differences for the final set (second iteration) generally lower than about ± 20%. By comparing the final (second iteration) mean spectra obtained using the TWM ( Figure 11E) and the SIA ( Figure 12A), or the corresponding percentage relative difference ( Figures 11F for the TWM and 12B for the SIA) it can be observed that the two proposed approaches lead to comparable results if the larger value of the time window is considered in the TWM.
The amplitude and the energy and frequency content of the final sets (second iteration) of the artificial accelerograms and of the target motion are also compared in Figures 13 (TWM) and 14 (SIA) where the total intensity I 0 and the mean period T m of each of the generated artificial accelerogram are plotted versus the corresponding peak acceleration PGA. In these figures the data plotted as a red rhombus represent the values of PGA, I 0 and T m of the selected motion (hereafter referred as target PGA, target I 0 and target T m ); the other data (grey dots) represent the N s pairs (PGA-I 0 or PGA-T m ) corresponding to each accelerogram of the artificial final set (second iteration) obtained using the TWM (Figure 13), with three different values of the time window (Δ̃= 0.25 s, Δ̃= 2 s, Δ̃= 4 s), and using the SIA (Figure 14). In the plots of Figures 13 and 14, the thin red lines detect the target PGA, I 0 and T m ; the thin black lines identify the corresponding F I G U R E 9 Cumulative zero-level up-crossings functions of the target motion and of the initial set of artificial accelerograms generated using the TWM with different time windows. mean values computed for each set of N s artificial accelerograms. The geometric distances between these lines quantifies the relative difference between the target and the mean values of PGA, I 0 and T m computed for each set of artificial accelerograms.
Specifically, with reference to the TWM (Figure 13), depending on the time windows, relative differences ranging between −19% and −25%, −31% and +4%, −10% and −18% can be observed between the mean value computed for the set of N s = 100 artificial accelerograms and the target PGA, I 0 and T m respectively; for the SIA ( Figure 14) these differences are equal to about −18%, −30% and −10%. Differently from the data previously commented with reference to the spectrum compatibility, Figure 13 clearly shows that the choice of the time window to be used in the TWM represents a crucial step. Specifically, the larger is the time window the larger is the difference between the target total intensity and the F I G U R E 1 0 Cumulative zero-level up-crossings function of the target motion and of the initial set of artificial accelerograms generated using the SIA.
corresponding average values of each set of accelerograms; conversely, the values of the time window do not significantly affect the relative differences in terms of PGA and T m .
Finally, as already observed with reference to the spectrum compatibility, the results obtained using the SIA and the TWM, with the larger values of the time window, are comparable. However, in both cases the procedure aimed to obtain the spectrum compatibility leads to final sets (second iteration) of accelerograms characterized by average total intensity, peak acceleration, and mean period farther from the target ones in comparison with the initial (0 th iteration) set. On the whole, the comparison of the plots given in the Figures 11-14 allows concluding that the better compromise in terms of spectrum, energy and frequency compatibility between the target motion and the sets of generated artificial accelerograms can be achieved by using the proposed TWM with the smallest value of the time window Δ̃= 0.25 s. Among the N s generated accelerograms those characterized by the lower differences with the target one in terms of PGA, I 0 and T m values could be detected in order to obtain an optimized sub-set of artificial accelerograms characterized by a very good compatibility with the target motion in terms of amplitude, energy and frequency content.
In this vein a sub-set of seven accelerograms with values of PGA and I 0 closest to the corresponding target ones have been selected among the final sets of the artificial motions generated using the TWM and the SIA; these subsets are represented with blue dots in Figures 13 and 14 A comparison between the target and the functions of the sets of N s = 100 accelerograms and the sub-sets of s = 7 artificial accelerograms obtained using the TWM (with Δ̃= 0.25 s) and the SIA are presented in Figure 15 in terms of Fourier amplitude spectra (Figures 15 A,B) and elastic response spectra assuming a 5% damping ratio (Figure 15 C,D). The blue lines in the plots of Figure 15 represent the mean functions computed for the sub-set of the selected seven accelerograms, while the black lines refer to the mean functions obtained considering the whole set of N s = 100 accelerograms.

Result of 1D non-linear site response analyses
The seismic site response analyses were carried out for the selected soil profile using the final sets (2 nd iteration) of N S = 100 generated artificial accelerograms as input motions. Some of the obtained results are presented in Figure 16 (TWM) and 17 (SIA) in terms of profiles of the maximum acceleration a max , stratigraphic amplification factor S S (computed as the ratio between a max and the peak acceleration of the rock outcropping motion) and peak shear strain γ max . Specifically, Figure 16 shows the results obtained using the N S input motions generated using the TWM with the three values of the time window already considered in Figure 13 (Δ̃= 0.25 s, 2 s, 4 s); similarly, Figure 17 refers to the results obtained using the N S input motions generated through the SIA. In these figures the grey area represents the envelope of the results whose mean values are represented by the black continuous line; in each plot, the dashed black lines define the confidence intervals of the obtained results computed, again, as the mean value plus/minus one standard deviation of the data distribution. For comparison, the results of the site response analyses obtained using the target accelerogram as input motion is also represented in Figures 16 and 17 with a red continuous line.
F I G U R E 1 1 Target elastic response spectrum (red solid line) together with the mean elastic spectra generated using the TWM, and corresponding percentage relative difference, for the sets obtained after the: 0 th (A,B), 1 st (C,D) and 2 nd (E,F) iterations required to achieve the spectrum compatibility.

F I G U R E 1 2
Comparison between the response spectrum of the target motion (red solid line) and the mean elastic spectra computed for the sets obtained after the 0 th , 1 st and 2 nd iterations required to achieve the spectrum compatibility by using the SIA (A); corresponding percentage relative difference (B).

F I G U R E 1 3
Values of the total Intensity I 0 (A-C) and of the mean period T m (D-F) versus the peak ground acceleration (PGA). Red lines for identifying the target values (red rhombus) and black lines for detecting the mean values of the set of N s = 100 artificial accelerograms (grey dots) generated using the TWM with different time windows; final set of seven accelerograms (blue dots).

F I G U R E 1 4
Values of the total Intensity I 0 (A-C) and of the mean period T m (D-F) versus the peak ground acceleration (PGA). Red lines for identifying the target values (red rhombus) and black lines for detecting the mean values computed for the set of N s = 100 artificial accelerograms (grey dots) generated using the SIA; final set composed by seven accelerograms (blue dots).
As a general remark, it can be noticed that the dispersion of the peak amplitude (PGA) and of the energy and frequency content of the artificial input motions (observed in the values of I 0 and T m already discussed with reference to Figures 13  and 14) is reflected in the large envelopes obtained for all the profiles represented in Figures 16 and 17.
The large variability of the results can be ascribed to two main combined effects: the possible coupling between the frequency content of the input motion and the vibration frequencies of the considered soil profile; the influence of soilnonlinearity which is related to the peak amplitude of the input motion and could affect the coupling phenomena.
However, it can be observed that all the confidence intervals (mean ± one standard deviation) represented in the plots of Figures 16 and 17 actually reveal a reduced variability of the computed responses which appear also unaffected by the values of the time interval adopted in the TWM and by the method (TWM or SIA) adopted for generating the artificial input motions.
At any depth (z), the responses computed for the target accelerogram are in a fair agreement with the average responses (and with the corresponding confidence intervals) predicted using all the sets of artificial accelerograms. A better agreement can be observed in the profiles relative to the peak acceleration and peak shear strain γ max while, the large dispersion of the PGA of the input motions justifies the larger differences in the profiles relative to the stratigraphic amplification factor S S (which do not represent a direct result of the site response analysis).
The results of the site response analysis obtained at the ground surface (z = 0) are presented in Figure 18 where the computed peak acceleration and stratigraphic amplification factor are plotted versus the PGA of the N s artificial input motions (grey dots); again, the red rhombus represents the data computed using the target accelerogram as input motion, while the blue dots represent the values of the best sub-set composed by seven accelerograms.
The results of the site response analysis obtained using only the sub-set of seven accelerograms selected among those generated using the TWM (with the lowest values of the time window) are presented in the plots of Figure 19. In these plots the light-blue shadowed area represents the envelope of all the computed values of a max , S S and γ max having the continuous blue line as mean profile. For comparison, the profile computed using the target accelerogram (red line) and the mean of the profiles computed using the whole set (N s ) of artificial accelerogram (black line) as input motion are also represented. As it can be observed, in all the cases the envelope of the data (blue area) defines narrow intervals of variation of a max , S S and γ max at any depth and the average profiles of each of these parameters (blue line) are always in a very good agreement with those computed for the target motion (red line). Also, a general better agreement is apparent in comparison with the mean profiles computed for the whole set of artificial accelerogram.

CONCLUDING REMARKS
In this paper, a novel procedure for generating fully non-stationary energy compatible artificial accelerograms has been presented. The proposed approach makes use of the evolutionary power spectral density (EPSD) function useful to generate artificial accelerograms having both the frequency and energy content very similar to that of a target accelerogram. The frequency content of the target accelerogram is first identified by analysing the variation over time of the number of zero-level up-crossings. Then, a set of stationary zero-mean filtered Gaussian subprocesses, with unitary variance, is introduced. In particular, a Gaussian subprocess is determined for each significant time interval in which there is a significative change in the number of zero level crossings. The power spectral density (PSD) function of each stationary subprocess is determined once the number of zero-level up-crossings and peaks is evaluated in the considered time interval. It has been also observed that the variation over time of the frequency content is significantly influenced by the geotechnical F I G U R E 1 6 Results of the 1D site response analysis carried out using the final set of the N s input motions generated using the TWM with different time windows: profiles of peak acceleration (A, D, G), stratigraphic amplification factor (B, E, H) and peak shear strain (C, F, I).
properties of the soil affecting the seismic response of the site of the recording station. Furthermore, the trend of the frequency variation is almost linear for rigid sites regardless the depth of the hypocentre. In this case only one PSD function could be determined.
Since the PSDs of the subprocesses are assumed with unitary variance, it is easy to calibrate the modulating function appropriately to guarantee a correspondence between the energy content of the artificial accelerograms and the target one.
In this paper, two methods, for the characterization of the modulating function, have been introduced, they are referred to as: the time window method (TWM) and the simplified intensity compatible approach (SIA). In particular, once the time variation of the cumulative intensity of the selected accelerogram is determined, the TWM captures the temporal evolution of the absolute amplitude by applying the central finite difference method. Conversely, the SIA provides a closed form expression for the modulating function which depends on the Arias intensity and on the strong motion duration of the F I G U R E 1 7 Results of the 1D site response analysis carried out using the final set of the N s input motions generated using the SIA: profiles of peak acceleration (A), stratigraphic amplification factor (B) and peak shear strain (C).

F I G U R E 1 8
Results of the 1D site response analysis carried out using the N s input motions generated through the TWM (2 nd iteration, Δ̃= 0.25 s) and the sub-set of seven motions (blue dots) having PGA and I 0 closest to the target values: peak acceleration (A) and stratigraphic amplification factor (B) at the ground surface (z = 0) versus the peak acceleration PGA of the input motion. Red lines for identifying the target values (red rhombus) and black lines for detecting the mean values of the set of N s = 100 artificial accelerograms (grey dots).

F I G U R E 1 9
Results of the 1D site response analysis carried out using the N s = 100 input motions generated through the TWM (2 nd iteration, Δ̃= 0.25 s) and the sub-set of seven motions having PGA and I 0 closest to the target values: profiles of peak acceleration (A), stratigraphic amplification factor (B) and peak shear strain (C). target motion. The main advantage in using the SIA is that the parameters of the modulating function can be evaluated by using empirical attenuation relationships, already available in the literature, instead of an acceleration time history.
The accuracy of the proposed procedure in generating sets of accelerograms with characteristics compatible with those of the target motion is illustrated in the paper through an application to equivalent-linear one-dimensional site response analyses, and by using both the proposed TWM and SIA. The artificial accelerograms are generated starting from a rock outcropping motion. The results of the site response analyses are presented and discussed in terms of profiles of peak accelerations, earthquake-induced shear strains and in terms of characteristics of the motion predicted at the ground surface. The obtained results pointed out that the proposed procedure is able to generate artificial accelerograms that catch both the energy variation and frequency content of the target motion and that, through a proper selection of a sub-set of the generated artificial motions, the prediction of the site seismic response is characterized by a reduced variability.
The characteristics of the generated artificial accelerograms pointed out that the TWM generally provides the better compromise in terms of amplitude, energy, frequency and spectrum compatibility with the target motion which should be defined through an acceleration time-history; conversely, the use of the proposed SIA, despite less accurate in terms of compatibility with the target accelerogram, requires the knowledge of few data, namely: Arias intensity, time extremes of strong motion duration and overall duration of the target accelerogram, as well as its number of zero level up-crossings and peaks. The first three parameters are used to calibrate the modulating function, the last two to define the PSD functions of the zero-mean stationary Gaussian stochastic subprocesses.