Pierre Gehl
email: p.gehl@brgm.fr

Darius M Seyedi
email: darius.seyedi@andra.fr

John Douglas
email: j.douglas@brgm.fr

Vector-valued fragility functions for seismic risk evaluation

Keywords: unreinforced masonry, uncertainties, fragility functions, ground-motion parameters, intensity measures, earthquake risk evaluation

This article presents a method for the development of vector-valued fragility functions, which are a function of more than one intensity measure (IM, also known as ground-motion parameters) for use within seismic risk evaluation of buildings. As an example, a simple unreinforced masonry structure is modelled using state-of-the-art software and hundreds of nonlinear time-history analyses are conducted to compute the response of this structure to earthquake loading. Dozens of different IMs (e.g. peak ground acceleration and velocity, response spectral accelerations at various periods, Arias intensity and various duration and number of cycle measures) are considered to characterize the earthquake shaking. It is demonstrated through various statistical techniques (including Receiver Operating Characteristic analysis) that the use of more than one IM leads to a better prediction of the damage state of the building than just a single IM, which is the current practice. In addition, it is shown that the assumption of the lognormal distribution for the derivation of fragility functions leads to more robust functions than logistic, log-logistic or kernel regression. Finally, actual fragility surfaces using two pairs of IMs (one pair are uncorrelated while the other are correlated) are derived and compared to scalar-based fragility curves using only a single IM and a significant reduction in the uncertainty of the predicted damage level is observed. This type of fragility surface would be a key component of future risk

Introduction

Seismic fragility analysis is inherently probabilistic as each of its constituents (e.g. hazard level, material properties, structural model and limit-state definition) is uncertain. The uncertainties can be split into two categories: aleatory (variability), which are inherently random, and epistemic, which are consequences of lack of knowledge. Most engineering applications are affected by both types of uncertainties. Although some researchers suggest that a clear distinction must be made between them (e.g. [START_REF] Paté-Cornell | Uncertainties in risk analysis: Six levels of treatment[END_REF], it is often difficult to determine whether an uncertainty is aleatory or epistemic from a modelling point of view [START_REF] Der Kiureghian | Aleatory or epistemic. Does it matter?[END_REF]. Practically speaking, this distinction is rather a choice of the modeller. Separating the uncertainties into these two categories helps us to determine which ones can be reduced [START_REF] Der Kiureghian | Aleatory or epistemic. Does it matter?[END_REF].

Two main sources of uncertainty, namely the variability in the ground motion and mechanical properties of the target building are considered in the following. A literature review shows several attempts at taking into account this random character on seismic response (e.g. [START_REF] Kwon | The effect of material and ground motion uncertainty on the seismic vulnerability curves of RC structure[END_REF][START_REF] Ellingwood | Quantifying and communicating uncertainty in seismic risk assessment[END_REF][START_REF] Fragiadakis | Fast performance uncertainty estimation via pushover and approximate IDA[END_REF]. [START_REF] Kwon | The effect of material and ground motion uncertainty on the seismic vulnerability curves of RC structure[END_REF] investigated the effects of strong-motion variability and random structural parameters on vulnerability curves of a reinforced concrete (RC) building and they conclude that the effect of strong-motion variability is much more important than the randomness in material parameters. [START_REF] Fragiadakis | Fast performance uncertainty estimation via pushover and approximate IDA[END_REF] introduced a static pushover-based method to estimate the uncertainty in seismic performance of a nine-storey steel-frame structure due to the variability of material properties.

Masonry constructions (the focus here) exhibit more variability in their mechanical properties compared to RC and steel frame structures. It is thus expected that such a ranking obtained for RC or steel structures cannot be applied directly to masonry constructions. This belief is supported by two recent studies that seek to assess the relative importance of the various sources of uncertainty in vulnerability asssessment for masonry buildings: seismic demand, structural capacity and limit-state definition [START_REF] Cattari | Probabilistic seismic damage scenario by mechanical models: The case study of Sulmona (Italy)[END_REF][START_REF] Pagnini | A mechanical model for the seismic vulnerability assessment of old masonry buildings. Earthquakes and Structures 2(1):25-Park Y, Ang A[END_REF]. They find, based on the capacity-spectrum approach, that variability in the structural capacity is the main contributor to the global dispersion of the fragility curves.

The estimate of the uncertainty accompanying calculated fragility can be integrated in a risk assessment analysis. As an example, [START_REF] Ellingwood | Quantifying and communicating uncertainty in seismic risk assessment[END_REF] proposed such a coupled risk assessment analysis. Estimates of the annual probability of exceeding pre-defined performance levels can then be computed.

On the other hand, a considerable source of epistemic uncertainty in seismic risk analysis is present in the construction of fragility curves. A fragility curve, by definition, represents the conditional probability of reaching a given damage state for a given hazard level characterized by one intensity measure (IM, also known as a ground-motion parameter).

The standard method to develop fragility curves neglects the uncertainty in the estimated damage caused by the use of a single IM, which cannot represent all characteristics of a ground motion. Recently, some efforts have been made to model the effect of several IMs on structural damage. [START_REF] Baker | A vector-valued ground motion intensity measures consisting of spectral acceleration and epsilon[END_REF] proposed a vector-based IM made of two parameters: spectral acceleration (SA) and epsilon, which is defined as the difference between the logarithms of observed SA and predicted SA at a given period. They showed that considering two IMs instead of one can improve the collapse probability calculation for a multi-degree-of-freedom structure. [START_REF] Kafali | Seismic fragility analysis: Application to simple linear and nonlinear systems[END_REF] used an alternative IM expressed by two parameters: earthquake magnitude and source-to-site distance. [START_REF] Rajeev | Increased accuracy of vector-IM-based seismic risk assessment?[END_REF] used SA at the first (T 1 ) and second eigenperiods (T 2 ) of the structure as the IMs. With reference to an example RC-frame structure, the accuracy of prediction of the seismic risk using the considered vector IM versus a conventional scalar IM is presented. They showed that an effective choice of T 2 leads to a better estimate of the seismic risk than that obtained employing a scalar IM consisting of SA at T 1 only, while reducing the associated dispersion in the estimate. However, for the studied example structure, the reduction is negligible in light of the effort required in switching from a scalar to a vector IM. [START_REF] Seyedi | Development of seismic fragility surfaces for reinforced concrete buildings by means of nonlinear timehistory analysis[END_REF] went the extra step in developing fragility functions for various damage states explicitly involving more than one IM (i.e. fragility surfaces) for an eight-storey RC building using nonlinear time-history analysis. [START_REF] Koutsourelakis | Assessing structural vulnerability against earthquakes using multidimensional fragility surfaces: A bayesian framework[END_REF] introduced a Bayesian framework to derive vector-valued fragility functions from the limited data available. He used a statistical-learning model based on logistic regression to assess the relative importance of different IMs: peak ground acceleration (PGA), root-mean-squared (RMS) intensity, [START_REF] Arias | A measure of earthquake intensity[END_REF] intensity and the power of the excitation spectrum at a period close to that of the structure.

The main goal of the present article is to introduce a complete framework to develop vector-valued fragility functions characterizing the physical vulnerability of buildings against earthquakes. A stochastic approach based on the results of several hundred nonlinear timehistory analysis is proposed. The selected IMs should be strongly correlated to the seismic response of the target structure. To this end, the damage level of an unreinforced masonry (URM) building is evaluated by the use of nonlinear numerical calculations. The whole structure is modeled by specific macro-elements representing the seismic behaviour of masonry panels. The maximum transient displacement at the top of the structure is used as the damage indicator. A statistical analysis is proposed to evaluate the correlation between different IMs and their correlation with the structural response. For each structural typology, a specific vector-valued IM can thus be defined. A method for deriving fragility surfaces is then proposed for the studied structure. Note that the derived surfaces are not proposed for use in practice because the example structure is a small box-type building without the geometrical complexity of real buildings. However, the presented method and discussion provide new insights on risk analysis of masonry structures.

General framework

The method described in this paper mainly relies on the exploitation of numerous nonlinear dynamic analyses because their use, as opposed to pushover curves and response spectra, allows study of the influence of many IMs. At first, a structural model of the studied building is considered and characterized. Modal analysis is performed to identify dynamic properties of the model. Damage limit states are determined through pushover analysis. Several variants of the structural model are also considered in order to account for variability in its mechanical properties. Then the model is subjected to many accelerograms that have been selected based on the seismotectonic context of the region of interest.

The results enable the study of the effects of many ground-motion characteristics on the building response, using various data mining techniques, such as: the variable clustering method, comparison of standard deviations and ROC (Receiver Operating Characteristics) analysis. The performance of single IMs can then be compared to combinations of two (or more) IMs, with respect to their ability to reduce the dispersion in fragility func-tions. Fragility curves, based on a lognormal distribution and a hybrid parameter, can finally be represented as fragility functions (or surfaces) with respect to two IMs. The proposed method is applied here, as an example, to a two-storey URM building. Each step is detailed in the following sections.

The proposed method aims to progress beyond the state of the art in the following ways.

The selection of adequate couples of IMs is based on a rigorous procedure and not solely on linear correlations (e.g. use of ROC analysis and creation of clusters of variables). The generally strong correlation between IMs had not been addressed in previous works (e.g. [START_REF] Seyedi | Development of seismic fragility surfaces for reinforced concrete buildings by means of nonlinear timehistory analysis[END_REF]. This issue is tackled thanks to the introduction of confidence bounds.

Besides, a review of the various probabilistic models that were recently proposed is made at the end of the paper, and the efficiency of different models is examined by considering the damage estimates obtained from the dynamic analyses. Finally, a set of vector-valued IMs are introduced. The use of a combination of these IMs enables the development of robust and analytical fragility functions that take the usual form of fragility curves and can be incorporated directly into existing seismic risk evaluation procedures.

Test structure

A numerical model was developed based on the results of an experiment at the University of Pavia [START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF] on a real-scale building. The exact geometry of the real building is considered, while the mechanical properties displayed in Table 1 were either taken from [START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF] (i.e. shear and compressive strength) or selected from common features of brick masonry and calibrated using the experimental pushover curve. The structure is 6.4 m high, with plan dimensions of 6 × 4.4 m. The mock-up consists in a single-cell building, without any internal walls. As axial forces play a crucial role in the strength of the masonry, concrete blocks were added on each floor to account for live loads: this is translated into the model by a vertical load of around 250 kN on each floor. The numerical model was built using the TREMURI code [START_REF] Lagomarsino | TREMURI program: Seismic analysis program for 3D masonry buildings[END_REF], which allows a masonry wall to be discretized into several components (vertical elements, named 'piers', horizontal ones, named 'spandrels', and rigid zones), through an equivalent-frame approach (see Figure 1). The behaviour of the masonry panel is represented by a macroelement developed by [START_REF] Gambarotta | On dynamic response of masonry panels[END_REF], and later improved by [START_REF] Penna | A macro-element procedure for the non-linear dynamic analysis of masonry buildings[END_REF] and [START_REF] Galasco | On the use of pushover analysis for existing masonry buildings[END_REF]. This nonlinear macromodel takes into account the common in-plane damage mechanisms of a masonry panel, namely: shear failure (sliding or diagonal cracking) and bending failure (tensile and compressive cracking). The out-of-plane failure mechanism of wall is not taken into account. Even though out-of-plane behaviour may play an important role in local damage, this mechanism is commonly neglected when the overall behaviour of a structure is being considered [START_REF] Calderini | In-plane strength of unreinforced masonry piers[END_REF].

[Fig. 1

about here.]

[Table 1 about here.] The building tested by [START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF] was subjected to a cyclic lateral loading on each of the two fac ¸ade walls: the fac ¸ade with the door (referred as the 'door wall', see Figure 1) is disconnected from the transverse walls in the experiment, thus allowing focus on the in-plane mechanism of a masonry panel. It was decided to model the same fac ¸ade wall to calibrate the model. Therefore, pushover analysis on the 'door wall' was compared to the results of the experiment and a good agreement, in terms of initial stiffness, maximum shear strength and yield displacements, between the model and the behaviour of the tested building was found (see Figure 2a). The experiment carried out by [START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF] consisted of a series of cyclic loadings of increasing amplitude: therefore only the extrema of each cyclic loading have been used to compare to the monotonic pushover analysis performed on the TREMURI model. Finally, the actual pushover curve for the whole building has also been computed using the results from the single wall analysis (Figure 2b). A force distribution according to the first modal shape has been used, whereas the pushover analysis on the fac ¸ade wall is based on a uniform force distribution (i.e. in the same conditions as the experiment).

[Fig. 2

about here.]

Modal analysis of the 3D model gives the first mode along the fac ¸ade of the building (referred to as the X-direction), with a period T = 0.149 s. Higher modes along the X-axis (torsion and opposite-floor displacements) are observed with periods around 0.05 s.

Pushover analysis of the building in the X-direction enables estimation of the ultimate displacement (d u = 18.8 mm) as well as the elastic-limit point (d y = 5.6 mm for a shear force of 210 kN), using an idealized bilinear curve that has been fitted to the actual pushover curve (see Figure 2b). The bilinear curve is calculated based on the Eurocode-8 criteria for the computation of ULS (ultimate limit-state) and DLS (damage limit-state).

The evaluation of damage states for the complete building is not trivial, as the criteria used to define the local damage level of a macroelement cannot be simply extrapolated to the global level. Numerous damage indicators are available in the literature, such as indices based on hysteretic-energy dissipation (e.g. Park and Ang, 1985), the determinant of the stiffness matrix or the maximum strain at various locations and floor displacements. The maximum transient displacement at the top of the structure is used in the following because its evaluation is straightforward and many correlations between floor displacement and damage states are available. [START_REF] Milutinovic | Risk-UE: An advanced approach to earthquake risk scenarios with applications to different European towns[END_REF] propose relations to link top displacements d y and d u to the EMS-98 (Council of Europe, 1998) damage scale (see Table 2). The studied structure has a predominant first mode, which corresponds to a uniform distribution of deformations along the building height. The pushover analysis has been carried out with an adaptive scheme based on the first modal shape, thus ensuring that the damage thresholds estimated from the static analysis and based on top displacement are adequate proxies for the floor drifts and a relevant indicator for both local and global damage. Still, it is noted that the choice of a given damage indicator may affect the results (e.g. identification of optimal IMs), as displacement-based damage measures are strongly correlated with 'peak' IMs (e.g. PGA and PGV), whereas energy-based indexes are more influenced by IMs representing the energy content of the signal (e.g. Arias intensity or cyclic damage parameters).

[Table 2

about here.]

Once the reference building model has been developed, the uncertainties in the mechanical properties of the brick masonry can be taken into account by deriving a set of variants with respect to the basic model. These variants are developed by randomly sampling material properties, which are assumed to follow a normal distribution. These distributions are defined by using the mean values in Table 1 and by setting coefficients of variation between 10 and 20% [START_REF] Rota | Processing Italian damage data to derive typological fragility curves[END_REF][START_REF] Rota | A methodology for deriving analytical fragility curves for masonry buildings based on stochastic nonlinear analyses[END_REF]. For each parameter, bounds (set at ±2σ ) prevent generation of unrealistic buildings. Twenty variant models are generated out of the parametric distribution. A Latin hypercube sampling scheme is used to capture the possible combinations with a reduced number of models. Thus, with 20 models, intervals of 5% can be defined and, for each parameter, an aleatory interval picked to sample the value for a given model. This procedure results in a total of 21 structures (one basic model plus 20 variants):

for each of them, a pushover analysis is performed to identify the boundaries of the EMS-98 damage states, which are specific to each variant.

Nonlinear dynamic analyses

Dynamic analyses of the models require a large set of strong-motion records to build robust fragility functions. A dataset of 777 accelerograms has been assembled, using both a sample of records from the European Strong-Motion Database [START_REF] Ambraseys | Dissemination of European Strong-Motion Data[END_REF] and synthetic signals generated using the non-stationary stochastic procedure proposed by [START_REF] Pousse | Non stationary stochastic simulation of strong ground motion time histories including natural variability: Application to the K-net Japanese database[END_REF]. The records correspond to shallow crustal earthquakes of magnitude (M w ) between 4.5 and 6.5 and epicentral distances between 10 and 100km . The signals also sample the five Eurocode-8 soil classes to account for local site variability. The consistency between real records and synthetic accelerograms has been checked for a wide selection of ground-motion parameters: the overall distribution of the two datasets (mean and standarddeviation) shows good agreement, except for duration-related IMs, which have a tendency to be underestimated in the synthetic signals used here [START_REF] Pousse | Non stationary stochastic simulation of strong ground motion time histories including natural variability: Application to the K-net Japanese database[END_REF].

The selected accelerograms are then applied at the base of the structure, along the Xdirection. The number of simulations has been optimized and the dataset of 777 accelerograms divided into 37 groups of 21 records each. Thus, within each group, it is possible to randomly assign each accelerogram to one of the 21 models. It is worth noting that the 37 groups were selected by ranking all accelerograms with respect to PGA (IM selected a priori to represent the 'damage potential' of the signal). This ensures that each group contains comparable signals, and as a result it can be assumed that all 21 models are subjected to accelerograms with similar intensity levels.

Finally, all nonlinear dynamic analyses are carried out and the maximum transient displacement at the top of the building is used to identify the different damage states, based on Table 2. The results of the analyses show that 63% of the simulation runs did not damage the structure, while about 17% of them induced collapse. Intermediate levels (slight to very heavy damage) are clearly under-represented: only 2% of the simulations stopped in the 'slight damage' state, 5% in 'moderate', 6% in 'extensive' and 7% in 'very heavy'. As a result, it appears that EMS-98 intermediate damage states are too narrow. The five damage states are thus merged into only two. As proposed by recent studies (Crowley et al, 2011), damage states from 'slight' to 'extensive' are considered as 'yield'; the damage states 'very heavy' and 'collapse' are also merged into one (near-collapse/collapse), as the numerical simulations near the ultimate deformation limits of the macroelements are usually less accurate and there is no benefit in setting a clear boundary between these damage states.

Studied IMs

For each of the accelerograms used in the simulations, fifty IMs are calculated, ranging from PGA to cyclic and duration parameters and elastic response spectral ordinates. One useful approach to present these parameters is to rank them according to their cross-correlation. To do this a variable clustering technique is performed using the TANAGRA software [START_REF] Rakotomalala | Tanagra: Un logiciel gratuit pour l'enseignement et la recherche[END_REF]. It consists of a succession of principal component analyses and the merging of the parameters with strong correlation into distinct class-variables, called 'clusters'. This procedure is helpful in reducing the many IMs to only three or four less-correlated parameters that are representative of the range of studied IMs (see Table 3). For each accelerogram, the ratios PGA/PGV and PGA/SA(0.15 s) are also computed, as such ratios roughly characterize the frequency content of the signals.

[Table 3 about here.]

Receiver Operating Characteristic (ROC) analysis

The efficiency of each of the IMs for the prediction of the damage state (none, yielding or collapse) of the structure must be evaluated. A convenient way of measuring the accuracy of a 'predictor' is to perform a ROC analysis on the data (e.g. [START_REF] Fawcett | ROC graphs: Notes and practical considerations for data mining researchers[END_REF]. This approach has previously been applied to earthquake risk evaluation by [START_REF] Koutsourelakis | Assessing structural vulnerability against earthquakes using multidimensional fragility surfaces: A bayesian framework[END_REF]. For a given IM, a set of data for 'model training' (e.g. development of a fragility model based on the chosen IM, which will then be compared to a second dataset in the 'validation' phase) is used. The first data set is used to build a fragility curve for a given damage state, based on a lognormal distribution, using the procedure described by [START_REF] Shinozuka | Statistical analysis of fragility curves[END_REF]. This approach consists of fitting the damage probability to a cumulative lognormal density function, the curve's parameters being estimated through a maximum-likelihood method, as the outcome of a Bernoulli experiment. Then, using a second set of (validation) data, a confusion matrix is built for different thresholds of damage probabilities. First, for each data point x i , the damage probability (for instance, yielding or collapse) is evaluated using the previously defined lognormal distribution: p i = P(d i ≥ DS|x i ). Then, for a range of probability thresholds p 0 varying from 0 to 1, p 0 and p i are compared: if p i ≥ p 0 , then the predicted damage d i is 1, and -1 otherwise. The predicted d i are compared to their actual value to assess whether the damage state is accurately estimated or not. Thus, for each value of the threshold p 0 , all data points can be classified into four possible combinations, which form the so-called confusion matrix: It can be observed that the most efficient IM will be the one that minimizes both 'misses' and 'false alarms'. Two measures that assess the quality of a prediction model are: sensitivity, which is the ratio of the number of hits to hits and misses, and specificity, which is the ratio of the number of correct rejections to correct rejections plus false alarms. The ROC curve is then constructed by plotting the different values of these ratios obtained for all thresholds p 0 , the vertical axis being the sensitivity and the horizontal one the complement of specificity (Figure 3).

-predicted d =
[Fig. 3

about here.]

One way to read these ROC curves is to check if they are close to the (0, 1) corner, which is the sign of a satisfactory model. In contrast, the 45 • line represents a perfectly random predictor and any ROC curve near that region implies a mediocre prediction model. As a result, whereas the shape of a curve cannot be easily exploited, a commonly-used indicator is the area under the ROC curve (AUC). AUC close to unity denotes an efficient model, whereas values close to 0.5 indicate little advantage over chance. However, it must be kept in mind that AUC is useful in evaluating the global performance of a parameter but two ROC curves may overlap locally. It has been proven that AUC is equivalent to the probability that this model will rank a randomly chosen positive instance (e.g. damage occurrence) higher than a randomly chosen negative instance (e.g no occurrence of damage) [START_REF] Fawcett | An introduction to ROC analysis[END_REF].

Performance of single IMs

AUC values are computed for each of the fifty IMs regressed over a lognormal distribution, and for each of the two considered damage states (yielding and collapse), to assess which ones are the most useful in predicting the effect of earthquake shaking for this structure.

Since no extra data are available to compose a validation set, a ten-fold cross-validation is performed, which is a commonly-used technique when data are limited. Two equal-sized bins are randomly formed out of the data (777 simulation results) and the probability model is built based on one bin, while the other is used to evaluate the model and perform the ROC analysis. This process is repeated ten times and the average value of AUC is used to obtain a performance indicator for each of the fifty IMs, for both damage states independently (see Table 4).

[Table 4

about here.]

A first comment to make is that the AUC are generally very high (close to unity), with only slight differences (around one decimal) between the IMs. The standard deviation of the lognormal distribution are also shown as a measure of the model dispersion, to demonstrate that the dispersion of the lognormal distribution (e.g. the slope of the fragility curve) does not totally determine the predictive power of the model. For example, Arias intensity leads to a fragility model with a fairly high standard deviation, even though its AUC is one of the highest.

It can be observed that one IM that is the best for both damage states cannot be found.

Acceleration-related parameters (e.g. PGA, A95, ARMS, ASI and cyclic damage parameters) and response spectral accelerations at short periods (up to 0.25 s) are efficient in predicting yielding, whereas, parameters like PGV, PGD, absolute durations and spectral accelerations at longer periods are more efficient for the prediction of collapse. It is also worth noting that the two groups of IMs that are emerging from the ROC analysis are more or less the same as those in clusters A and B identified by the variable clustering analysis performed above. This observation reinforces the idea that a combination of two IMs may improve the prediction of damage, as the first could be useful in predicting slight damage, and the second more efficient for higher damage states.

The spectral accelerations at the two periods of interest for the studied structure (0.15 s and 0.05 s) are useful for the prediction of damage, especially for yield, yet these parameters are outranked by SA(0.25 s) and SA(0.5 s), even for yield. This may be explained by the rapid loss of initial stiffness of the masonry model, even for light damage. This issue has also been raised by [START_REF] Bommer | The influence of strong-motion duration on the seismic response of masonry structures[END_REF], where the authors propose to use a spectral acceleration averaged over a period interval of width several times the initial eigenperiod.

Linear combinations of two IMs

In the previous section, ROC analysis showed that two IMs might lead to more accurate assessment of both yield and collapse. For slight damage, an appropriate IM (IM1) could be chosen from cluster A (see Figure 3), whereas an IM (IM2) selected from cluster B should accurately predict heavier damage. A hybrid IM is, therefore, proposed using a linear combination of IM1 and IM2:

log X i = α 1 α 1 + α 2 log IM1 i + α 2 α 1 + α 2 log IM2 i (1)
where α 1 and α 2 are regression coefficients. This new parameter can now be considered as a single variable and the probability of reaching or exceeding a given damage state given X can then be expressed via the usual lognormal distribution [START_REF] Shinozuka | Statistical analysis of fragility curves[END_REF]:

P(damage ≥ DS|X) = 1 2 1 + erf log X -µ σ √ 2 (2)
where µ is the mean, σ the standard deviation, DS damage state and erf is the error function.

Then, by expressing X as a function of IM1 and IM2 and by introducing

β 1 = α 1 (α 1 +α 2 )σ √ 2 , β 2 = α 2 (α 1 +α 2 )σ √ 2 and β 0 = µ σ √
2 , the following relation is obtained:

P(damage ≥ DS|IM1, IM2) = 1 2 [1 + erf(β 1 log IM1 + β 2 log IM2 -β 0 )] (3)
Finally, the coefficients β i are evaluated through the maximization of the following like- lihood function:

L = n ∏ i=1 P i y i (1 -P i ) 1-y i (4)
where y i is the so-called binomially distributed variable and equals unity when the damage state is reached, and null otherwise and P i is the probability of reaching or exceeding the damage state given (IM1 i , IM2 i ).

Using the relations presented above, several combinations of IMs are tested and these new fragility models are evaluated through a ROC analysis, using again a ten-fold crossvalidation scheme. The results obtained for some pairs of IMs are presented in Table 5.

[Table 5 about here.]

Comparing the AUC found for two IMs (Table 5) with those in Table 4 shows a general improvement in the prediction of damage states. Also, the standard deviations have been reduced for both damage states, in comparison to single IM models, thus showing the performance of multivariate fragility functions in terms of uncertainty reduction. It is now possible to accurately predict both damage states (yield and collapse) at the same time, by using a combination of IMs from clusters A and B. Nonetheless, it is not strictly correct to speak of 'hybrid' parameters (e.g. a function of two IMs) because the coefficients in the linear combination are specific to each damage state, depending on the relative influence of each parameter. It is thus preferable to refer to them as vector-valued parameters, whose scalar components can be linearly combined in a regression model.

Table 5 indicates that the optimal combination of IMs would be the couple SA(0.25 s)-SA(0.5 s). However, the spectral acceleration at the first vibration mode (i.e. 0.15 s) could also be used, as this choice seems more in line with the structural study. One also has to keep in mind that other considerations can influence the choice of the IM, such as its ability to be predicted via GMPEs (e.g. [START_REF] Douglas | Consistency of ground-motion predictions from the past four decades: Peak ground velocity and displacement, Arias intensity and relative significant duration[END_REF]. It is assumed that the very slight differences in AUC values displayed in Table 5 are not significant; therefore, any couples of IMs presented in this table can be considered for use, as they are at least as efficient as any single IM (see Table 4).

Finally, it is interesting to observe that some parameters that have little effect on the structural response when considered alone [e.g. the ratios PGA/PGV or PGA/SA(0.15 s)]

can have a positive impact on the prediction model when combined with a second IM. For example, the couple (PGV, PGA/PGV) has a greater AUC (Table 5) than just PGV or PGA alone (Table 4). [START_REF] Kwon | The effect of material and ground motion uncertainty on the seismic vulnerability curves of RC structure[END_REF] also stress the role of the PGA/PGV ratio as an indicator of the dynamic characteristics of ground-motion, since impulsive records with a low PGA/PGV ratio are usually characterized by high energy contents and low frequencies.

Development of bivariate fragility functions

Considering Table 5, there is a wide choice of adequate couples of IMs to use as variables in fragility functions. The following section details the development of fragility functions and their implications in terms of dispersion reduction, through two distinct examples that highlight the importance of accounting for the correlations between IMs.

Case of uncorrelated IMs

Since the couple (PGV, PGA/PGV) is efficient and these two parameters are almost orthogonal (correlation coefficient equal to -0.1548), they are chosen for the first example fragility surface (see Figure 4), using the formulation of Equation 3. The regression coefficients obtained through the likelihood method are given in Table 6.

[Fig. 4 about here.]

[Table 6 about here.]

The choice of two parameters that are uncorrelated or almost orthogonal is convenient, as it allows evaluation of the vulnerability on the whole 2D space defined by PGV and PGA/PGV, even in the corners that contain high values of one parameter and low values of the second. To study the effect of the second parameter, single-variable fragility curves are plotted by keeping the second IM constant and then compared with a fragility curve built by considering only the first parameter (see Figure 5).

[Fig. 5 about here.]

First, it can be noticed that the use of a second parameter (here the ratio PGA/PGV) induces a steady increase in the damage probability. For instance, in the case of small PGVs (PGV around 0.05 m/s), the single-variable curve shows almost zero probability of damage; however, the fragility surface indicates that if this signal is accompanied by high accelerations (e.g. high PGA/PGV ratios around 20 or 30), then there is actually a non-negligible probability of damage (around 0.25). A crucial observation is that the single-variable curve is flatter than those extracted from the surface: this confirms that fragility models built with a single IM have a greater dispersion than surfaces and they are associated with higher epistemic uncertainty because earthquake signals are only characterized by a single number.

Case of correlated IMs

The IMs chosen above are convenient to demonstrate the benefit of fragility surfaces because they are almost orthogonal. However, the ratio PGA/PGV is not a standard IM and its estimation by ground-motion prediction equations in the framework of seismic hazard assessment is not easy. Therefore, a fragility function with respect to SA at different periods is derived here. These IMs are more convenient for joint estimation and they are well adapted to vector-valued probabilistic seismic hazard assessment (VPSHA) proposed by [START_REF] Bazzurro | Vector-valued probabilistic seismic hazard analysis (VPSHA)[END_REF]. The SA at the first eigenperiod (0.15 s) and a second one at a longer period are chosen, based on the results from the ROC analysis (see Table 5). For this second period 0.5 s was chosen as it has the ability to represent the response of the building once its initial stiffness has been reduced due to damage. The AUC for collapse for SA(0.5 s) is the highest of all of the IMs considered (Table 4) hence it was preferred over SA at other periods. The optimum IMs for the construction of fragility curves depends on the structure.

For example, for a building with a higher eigenperiod (e.g. a taller building) SAs at longer periods will probably be more effective as a basis of fragility curves. Similarly the level of nonlinearity in the building response will also affect what periods are optimum.

The spectral accelerations at 0.15 and 0.5 s are, however, strongly correlated (R = 0.81) and caution must be taken in the construction of the surfaces and their interpretation. Due to the strong correlation between the two IMs, the data points do not cover the whole space defined by SA(0.15 s)-SA(0.5 s) and one can question the relevance of defining a fragility model for extreme values (e.g. low IM1 and high IM2). Using the data points, a regression analysis is performed between SA(0.15 s) and SA(0.5 s), allowing a median line representing the linear relation between the two parameters to be plotted, as well as 95% confidence intervals (see Figure 6). The median line is represented by the equation log SA(0.5) = log SA(0.15) -0.7236 ± δ , with the 5 -95% lower and upper bounds estimated with δ = 1.417. Iso-probability lines are also represented for both damage states and the grey area between the confidence intervals can be considered as the validity domain for the fragility surfaces. Like in the previous example, the fragility surfaces are built assuming a lognormal distribution (Equation 3); the corresponding regression parameters are listed in Table 6.

[Fig. 6 about here.]

Finally, as described previously, single-variable fragility curves can be plotted with respect to the first IM, by fixing the second one. However, in the present case, the two IMs are correlated and plotting the fragility function using one IM while keeping the other one constant may lead to erroneous conclusions, especially for extreme values. Thus, one solution may be to cut 'slices' along the affine lines within the validity domain defined on Figure 6, to account for the correlation between the two IMs. The results for yield damage state are shown on Figure 7, for both IMs, and for a few correlation hypotheses (median relation between the two IMs, and 95% upper and lower bounds). Again, the effect of the introduction of a second IM can be observed, which allows for a clearer identification of the damage states, as shown by the steeper slopes of the fragility functions extracted from the surfaces.

[Fig. 7 about here.]

Discussion

The fragility functions derived here result from several trials on which probability distribution to adopt for the surfaces. First, a logistic regression, as suggested by [START_REF] Koutsourelakis | Assessing structural vulnerability against earthquakes using multidimensional fragility surfaces: A bayesian framework[END_REF] was tested. This approach contains some drawbacks: the logistic function has a support on both positive and negative real numbers, which leads to a possibility of having nonzero probabilities even when IMs are equal to zero (this high kurtosis can be an problem when a risk assessment is performed for low hazard sites). It was attempted to overcome this problem by using the log-logistic distribution (i.e. logistic distribution of the logarithm of a variable), which constrained the probabilities at the origin. However, a ROC analysis of the models showed that the distribution presented in Equation 3 (based on a lognormal formulation) yielded slightly higher AUC values than the log-logistic distribution.

Gehl et al ( 2011) proposed an alternative solution, using a kernel density estimation procedure (or Parzen-Rosenblatt method), which allows a non-parametric definition of the probability density and is, thus, applicable to higher-dimension problems [START_REF] Parzen | On estimation of a probability density function and mode[END_REF].

The main merit of this approach is that the shape of the cumulative density function is not biased by the choice of the distribution and it can account for all the specificities of a bivariate probability density. However, the Parzen-Rosenblatt formulation relies strongly on a 'smoothing' parameter (e.g. a standard-deviation for each data point), which greatly influences the slopes of the surfaces and is also specific to the data used (quantity of points and distribution). Moreover, the kernel density estimation cannot be represented by a simple analytical relation and it is numerically estimated, which can be a problem if the fragility function is to be used in risk assessments.

[Fig. 8 about here.]

In Figure 8, for the same data [yield damage state based on the couple (PGV, PGA/PGV)] the iso-probability lines for the four models discussed in this section are plotted. The logistic distribution is not very adapted to the studied problem, as it fails to accurately represent the damage for low PGA/PGV (near the X-axis of the plot). The kernel density estimation approach has the ability to best represent the results, yet this empirical distribution has to be smoothed to reduce local irregularities, which has the non-desired effect of flattening the slope of the surface (e.g. widening of the 16-84% percentile zone). Both log-logistic and lognormal distribution are able to fit the general shape of the kernel density estimation curve and they are very similar. The lognormal distribution was finally selected to build the fragility surfaces because, as explained above, the ROC analysis of this model gave slightly better results than the log-logistic one. Also, this distribution is consistent with the lognormal distribution commonly used for single-IM fragility curves.

Conclusions

Current methods used to evaluate seismic risk are based on fragility curves representing the ground motion by a single IM (e.g. PGA). Different types of uncertainty affect the results of such evaluations. Three sources of aleatory and epistemic uncertainty are addressed in the present work. The variability of the ground motion and randomness of the material properties of the buildings are taken into account through a large number of nonlinear time-history calculations. Several hundred accelerograms and several variants of the target building are used to this end. Moreover, representing earthquake shaking by a single IM introduces epistemic uncertainties in the risk evaluation as such a model cannot fully represent the effect of an earthquake on the response of the structure. It is shown that an increase from one to two IMs leads to a significant reduction in the scatter in the fragility function and consequently a potential reduction in the uncertainty in evaluated earthquake risk.

Suggestions for future research are the following. It would be useful to investigate the effect of the out-of-plane failure on fragility curves of URM buildings. In large-scale seismic risk assessments, the behaviour of a typology of buildings is generally studied. The effect of geometrical variability on the vulnerability of a given typology must be taken into account by introducing a population of different buildings belonging to the same typology. To make use of the vector-valued fragility function like those proposed here within risk evaluations, it would be necessary to conduct a VPSHA [START_REF] Bazzurro | Vector-valued probabilistic seismic hazard analysis (VPSHA)[END_REF] so that the joint probability of surpassing thresholds for two (or more) IMs was correctly estimated. VPSHA has yet to become a standard procedure in risk evaluations but thanks to recent developments of the method (e.g. [START_REF] Bazzurro | Efficient approach to vector-valued probabilistic seismic hazard analysis of multiple correlated ground-motion parameters[END_REF] this approach is likely to become more widespread.

List of Figures

1 Numerical model of the building tested by [START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF] Fig. 6 Fragility surfaces using the couple [SA(0.15 s), SA(0.15 s)] as earthquake descriptors, for damage states 'yield' (left) and 'collapse' (right). The solid black line is the result of linear regression between SA(0.15 s) and SA(0.5 s), and the dotted lines are the 95% confidence intervals. The grey area is the proposed validity domain for the fragility surfaces. The blue dots are the data points used in the study. The colored curves are the iso-probability lines for probability levels 5, 16, 50, 84 and 95%. Fig. 8 Comparison of the four probability models considered here (logistic, log-logistic, kernel density estimation and lognormal). The blue dots represent the data that have not reached the yield damage state, and the red crosses the points that induced damage. The solid black line represents the 50% probability of reaching the damage state, and the dotted lines the 16 and 84% percentiles. Absolute significant duration ARMS SA(T = 0.5 to 5 s) ASI Cyclic damage parameters SA(T = 0.05 to T = 0.25 s) Note: The spectral accelerations (SA) were computed for a 5% damping ratio, for the following periods: 0.05 (second eigenperiod), 0.075, 0.1, 0.15 (first eigenperiod), 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2 and 5 s. A95 is the level of acceleration that contains 95% of the Arias intensity [START_REF] Sarma | An evaluation of strong motion records and a new parameter A95[END_REF]. SLOPE75 (and SLOPE95) are the slope of the Husid plot (e.g. cumulative Arias intensity over time) between 5% and 75% (and 95%) of the total Arias intensity. ARMS (root-mean-square acceleration) is the square-root of the integral of squared acceleration over time. ASI (Acceleration spectral intensity) is the integral of SA between two periods (here, 0.1 s and 0.5 s). Cyclic damage parameters are obtained through the sum of the squared amplitude of all half-cycles, which are obtained through three different counting methods (rainflow and peak counting, with or without non-zero crossings) [START_REF] Hancock | The effective number of cycles of earthquake ground motion[END_REF]. Nb of effective cycles is the same as the cyclic damage parameters, except that the half-cycles amplitudes are normalized by the amplitude of the largest half-cycles in the signal [START_REF] Hancock | The effective number of cycles of earthquake ground motion[END_REF]. Absolute bracketed duration is the length of interval between the first and last time the ground acceleration exceeds a threshold value (here 0.05 g). Absolute uniform duration is the total length of time for which ground acceleration exceeds a threshold value (here 0.05 g). Absolute significant (effective) duration is the length of interval between when Arias intensity first exceeds a threshold value (here 0.01 m/s) and the time when Arias intensity first exceeds total Arias intensity of record minus some threshold value (here 0.125 m/s) [START_REF] Bommer | The effective duration of earthquake strong motion[END_REF]. Relative bracketed duration is the length of interval between the first and last time the ground acceleration exceeds 5% of PGA. Relative uniform duration is the total length of time for which ground acceleration exceeds 5% of PGA. Relative significant duration is the length of interval between when Arias intensity first exceeds 5% of total value and time when Arias intensity first exceeds 95% of total value [START_REF] Trifunac | A study on the duration of strong earthquake ground motion[END_REF]. Note: Only the cyclic damage parameters based on rainflow counting are listed (peak counting methods yield almost identical results). Also, IMs that did not prove to be efficient for damage prediction (AUC values around 0.5): relative bracketed duration, relative uniform duration, (relative and absolute) significant duration, equivalent number of cycles, ratios PGA/SA(0.15 s) and PGA/PGV, are not listed.

Table 5 Computed AUC values for the most efficient couples of IMs for the two damage states (yield and collapse). In the fifth column, the linear correlation coefficient R between the two IMs is given. The standard deviations σ of the lognormal distribution for both damage states are also listed in the final two columns.

  1 and actual d = 1: true positive or 'hits', -predicted d = 1 and actual d = -1: false positive or 'false alarms', -predicted d = -1 and actual d = 1: false negative or 'misses', -predicted d = -1 and actual d = -1: true negative or 'correct rejections'.

Fig. 1

 1 Fig.1Numerical model of the building tested by[START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF].

Fig. 3 Fig. 4 Fig. 5

 345 Fig. 2 a) On the left, pushover results for a fac ¸ade wall and comparison between the TREMURI model and the experimental data by Magenes et al (1995). b) On the right, modal pushover on the whole building model and corresponding bilinear idealization.

Fig. 7

 7 Fig.7Single-variable fragility curve (solid black lines) compared to slices of a fragility surface (dotted lines), for both IMs.

  . . . . . . . 23 2 a) On the left, pushover results for a fac ¸ade wall and comparison between the TREMURI model and the experimental data by[START_REF] Magenes | Static testing of a full scale, two-story masonry building: test procedure and measured experimental response, in experimental and numerical investigation on a brick masonry building prototype[END_REF]. b) On the right, modal pushover on the whole building model and corresponding bilinear idealization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Example of ROC curves for three ground-motion parameters for the yield damage state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4 Example of a fragility surface for 'Collapse' state, using the couple (PGV, PGA/PGV) as earthquake descriptors. . . . . . . . . . . . . . . . . . . . . 26 5 Single-variable fragility curve (solid black line) compared to slices of a fragility surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6 Fragility surfaces using the couple [SA(0.15 s), SA(0.15 s)] as earthquake descriptors, for damage states 'yield' (left) and 'collapse' (right). The solid black line is the result of linear regression between SA(0.15 s) and SA(0.5 s), and the dotted lines are the 95% confidence intervals. The grey area is the proposed validity domain for the fragility surfaces. The blue dots are the data points used in the study. The colored curves are the iso-probability lines for probability levels 5, 16, 50, 84 and 95%. . . . . . . . . . . . . . . . . . . . 28 7 Single-variable fragility curve (solid black lines) compared to slices of a fragility surface (dotted lines), for both IMs. . . . . . . . . . . . . . . . . . 29 8 Comparison of the four probability models considered here (logistic, loglogistic, kernel density estimation and lognormal). The blue dots represent the data that have not reached the yield damage state, and the red crosses the points that induced damage. The solid black line represents the 50% probability of reaching the damage state, and the dotted lines the 16 and 84% percentiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Table 2

 2 Correlation between EMS-98 damage states and the transient top displacement, for URM buildings, according to[START_REF] Milutinovic | Risk-UE: An advanced approach to earthquake risk scenarios with applications to different European towns[END_REF].

	Damage state	Limit-displacement equation	Limit for d ( mm)
	Slight	d = 0.7d y	3.94
	Moderate Extensive Very heavy Collapse	d = 0.7d y + 0.05(0.9d u -0.7d y ) d = 0.7d y + 0.2(0.9d u -0.7d y ) d = 0.7d y + 0.5(0.9d u -0.7d y ) d = 0.9d u	4.59 6.54 10.43 16.92

Table 3

 3 Clusters generated from the set of accelerograms.

	Cluster A	Cluster B	Cluster C
	PGA	PGV	Relative bracketed duration
	Arias (1970) intensity	PGD	Relative uniform duration
	A95	Absolute bracketed duration	Relative significant duration
	SLOPE75	Absolute uniform duration	Nb of effective cycles
	SLOPE95

Table 4

 4 Computed AUC values for the IMs and for the two damage states (yield and collapse). In italics, parameters with AUC over 0.95, and in bold IMs with AUC above 0.98. The standard-deviation σ of the lognormal distribution for both damage states is also given.

	IM	AUC	AUC	σ	σ
		(Yield)	(Collapse)	(Yield)	(Collapse)
	PGA	0.9856	0.9680	0.2885	0.4957
	PGV	0.9674	0.9786	0.4211	0.3467
	PGD	0.8333	0.8804	1.4637	1.4694
	Arias intensity	0.9817	0.9834	0.6047	0.6846
	A95	0.9883	0.9710	0.2626	0.4943
	ARMS	0.9832	0.9677	0.3296	0.5174
	ASI	0.9771	0.9683	0.4291	0.5413
	SLOPE75	0.9878	0.9679	0.5826	1.1240
	SLOPE95	0.9870	0.9691	0.5887	1.0710
	SA(0.15 s)	0.9682	0.9536	0.4530	0.6184
	SA(5.0 s)	0.8576	0.8980	1.2480	1.0102
	SA(2.0 s)	0.9208	0.9492	0.8517	0.6567
	SA(1.75 s)	0.9274	0.9515	0.7995	0.6110
	SA(1.5 s)	0.9306	0.9623	0.7487	0.5491
	SA(1.25 s)	0.9348	0.9665	0.7136	0.4998
	SA(1.0 s)	0.9410	0.9728	0.6328	0.4490
	SA(0.75 s)	0.9610	0.9799	0.4966	0.3498
	SA(0.5 s)	0.9814	0.9876	0.3291	0.2585
	SA(0.25 s)	0.9873	0.9712	0.2724	0.4483
	SA(0.1 s)	0.9571	0.9353	0.5161	0.7917
	SA(0.075 s)	0.9633	0.9427	0.5108	0.7808
	SA(0.05 s)	0.9714	0.9496	0.4700	0.6981
	Duration bracketed absolute	0.9284	0.9264	0.8970	0.7662
	Duration uniform absolute	0.9746	0.9781	0.6607	0.4997
	Cycle damage paramater (rainflow)	0.9722	0.9641	0.8484	1.0554

Table 6

 6 Regression coefficients computed for the couples (PGV, PGA/PGV) and [SA(0.15 s), SA(0.5 s)] for damage states 'yield' and 'collapse'.

	IM1	IM2	AUC	AUC	Pearson R	σ	σ
			(Yield)	(Collapse)		(Yield)	(Collapse)
	SA(0.25 s)	SA(0.5 s)	0.9932	0.9916	0.8798	0.2007	0.2424
	PGA	SA(0.5 s)	0.9928	0.9902	0.8674	0.1927	0.2484
	ASI	SA(0.5 s)	0.9906	0.9910	0.8766	0.2554	0.2572
	SA(0.05 s)	SA(0.5 s)	0.9909	0.9903	0.7789	0.2308	0.2575
	SA(0.15 s)	SA(0.5 s)	0.9904	0.9891	0.8116	0.2293	0.2537
	AI	SA(0.5 s)	0.9874	0.9920	0.7731	0.3532	0.2966
	ARMS	SA(0.5 s)	0.9896	0.9898	0.7979	0.2291	0.2686
	Cyclic damage (rainflow)	SA(0.5 s)	0.9879	0.9906	0.7073	0.3139	0.2831
	PGV SA(0.25 s)	PGA/PGV PGV	0.9923 0.9916	0.9845 0.9848	-0.1548 0.8516	0.1313 0.2233	0.2234 0.3008
	SA(0.25 s)	SA(1.0 s)	0.9897	0.9864	0.8111	0.2340	0.3046
	PGA	SA(1.0 s)	0.9899	0.9854	0.8168	0.2409	0.3140
	PGA	PGV	0.9912	0.9840	0.8720	0.2088	0.3118
	ASI	PGV	0.9888	0.9862	0.7988	0.2495	0.3043
	PGA	PGA/PGV	0.9911	0.9834	0.1450	0.3538	0.7883
	AI	PGV	0.9860	0.9874	0.7830	0.4044	0.4168
	SA(0.25 s)	D uni. abs.	0.9898	0.9836	0.8064	0.2866	0.3702
	SA(0.05 s)	PGV	0.9896	0.9831	0.7897	0.2413	0.3268
	ARMS	SA(1.0 s)	0.9867	0.9858	0.7598	0.2742	0.3382
	PGA	D uni. abs.	0.9883	0.9833	0.7397	0.2994	0.3733
	ARMS	PGV	0.9888	0.9824	0.8257	0.2379	0.3284
	SA(0.15 s)	PGV	0.9873	0.9835	0.8126	0.2541	0.3178
	AI	SA(1.0 s)	0.9831	0.9872	0.7331	0.5399	0.4699
	ARMS	D uni. abs.	0.9867	0.9836	0.6846	0.3397	0.3858

Acknowledgements

The work presented in this article has been partially supported by the French Research National Agency (ANR) through the PGCU-2007 program (Project: A mechanical approach for evaluating seismic vulnerability of masonry structures, EVSIM, under grant number: ANR-07-PGCU-002), and the PERPETUATE (Performance-based approach to earthquake protection of cultural heritage in European and Mediterranean countries) project of the EC-Research Framework Programme FP7. In addition, we thank Guillaume Pousse, for his computer program to simulate ground motion, and Prof. Sergio Lagomarsino, for sharing with us the research and development version of TREMURI. Finally, we thank the two anonymous reviewers for their constructive comments on a previous version of this article.