7 Seismic fragility models provide a probabilistic relation between ground motion intensity 8 and damage, making them a crucial component of many regional risk assessments. Estimat- 9 ing such models from damage data gathered after past earthquakes is challenging because 10 of uncertainty in the ground motion intensity the structures were subjected to. Here, we 11 develop a Bayesian estimation procedure that performs joint inference over ground motion 12 intensity and fragility model parameters. When applied to simulated damage data, the pro- 13 posedmethodcanrecoverthedata-generatingfragilityfunctions, whilethetraditionallyused 14 method, employing fixed, best-estimate, intensity values, fails to do so. Analyses using syn- 15 thetic data with known properties show that the traditional method results in flatter fragility 16 functions that overestimate damage probabilities for low-intensity values and underestimate 17 probabilities for large values. Similar trends are observed when comparing both methods on 18 real damage data. The results suggest that neglecting ground motion uncertainty manifests 19 in apparent dispersion in the estimated fragility functions. This undesirable feature can be 20 mitigated through the proposed Bayesian procedure. 21 ∗Corresponding author: bodenmann (at) ibk.baug.ethz.ch 1 INTRODUCTION 22 Procedures for earthquake risk assessments often rely on probabilistic relations between 23 ground shaking intensity and physical damage to elements of the built environment. A 24 fragility function plots the probability of reaching or exceeding a certain damage limit state 25 for increasing ground motion intensity measure ( IM) levels. The fragility function parame- 26 ters are derived from expert opinion (e.g., Rojahn and Sharpe, 1985), or estimated from data 27 that is either generated through structural analyses of building archetypes, the so-called me- 28 chanical or analytical approach (e.g., Mosalam et al., 1997), or gathered after an earthquake 29 event, the so-called empirical approach (e.g., Basöz and Kiremidjian, 1998). 30 Besides challenges in quality and coverage of post-earthquake damage surveys, empiri- 31 cal fragility modeling also requires assumptions on the IMvalues that caused the observed 32 damage (Rossetto and Ioannou, 2018). These values are only known at locations of seismic 33 network stations, and are uncertain at other locations. Thus, if a group of damaged build- 34 ings is observed, we do not know a-priori whether they were damaged because they were 35 particularly vulnerable, because the ground shaking was particularly strong, or because of 36 both. This study aims to provide further insights into this causality dilemma and proposes 37 a new methodology to estimate fragility function parameters. 38 One common approach to address this causality dilemma assumes fixed, best-estimate, 39 IMvalues for parameter estimation (Rossetto et al., 2014). The fixed values represent best- 40 estimates that are typically chosen as the median of the IMdistribution conditional on 41 available data from seismic network stations (e.g., from a shake map). This fixed IMap- 42 proach served as the basis for the vast majority of available empirical fragility functions, 43 including early (e.g., Basöz and Kiremidjian, 1998) and recent studies (e.g., Rosti et al., 44 2021). By considering fixed IMvalues, such an approach does not account for potential 45 deviations of the actual IMfrom this best-estimate. The sensitivity of estimated fragility 46 2 function parameters with respect to IMuncertainty was examined by Ioannou et al. (2015) 47 through simulated case studies. They found that unrealistically dense seismic networks are 48 required to recover the “true” fragility parameters used to generate the data. For realistic 49 network configurations, however, the approach based on fixed IMvalues led to substantial 50 differences between estimated and “true” parameters. Approaches to account for IMuncer- 51 tainty vary in the existing literature. For example, King et al. (2005) focused on damage 52 survey data collected within a limited distance to a seismic network station. Lallemant et al. 53 (2015) suggested weighting the observations according to the inverse of the remaining IM 54 variance, thus assigning higher weight to data collected close to stations. Ioannou et al. 55 (2015) simulated numerous IMvalues at the survey locations and performed parameter esti- 56 mation for each simulation run. The resulting set of estimated parameters was then assumed 57 to reflect IMuncertainty. 58 While aiming to account for IMuncertainty, the above approaches miss a crucial compo- 59 nent of the causality dilemma, which is best explained by considering the inverse problem, 60 where the true fragility function parameters are available. In this hypothetical case, studied 61 by Pozzi and Wang (2018), the damage observations serve as noisy measurements of the 62 uncertain IMvalues. Consequently, these observations can be used to obtain a constrained 63 estimate of the IMvalues that caused the damage. 64 In cases of practical interest, both the fragility parameters and the IMvalues are uncer- 65 tain. This study presents a Bayesian approach to estimate their joint posterior distribution, 66 which can account for uncertainty in both the IMvalues and the fragility parameters. We 67 start with an overview of fragility function definitions, followed by mathematical descriptions 68 of the proposed Bayesian parameter estimation approach and the traditional, fixed IM, ap- 69 proach. Both methods are first explained on a simplified example and then compared using 70 data from the 2009 L’Aquila (Italy) earthquake. 71 3 FRAGILITY FUNCTIONS 72 Empirical fragility modeling aims to estimate the fragility function parameters using ob- 73 served or estimated values of IMs and building damage states. Damage states are obtained 74 from post-earthquake building inspections, where experts describe the building characteris- 75 tics and the sustained damage based on pre-specified guidelines and data-recording forms. 76 Information on the building characteristics pertain to the material and the type of the load 77 resisting system, as well as geometrical attributes such as the number of storeys. Build- 78 ings with similar characteristics are grouped into a ‘building class’ with a corresponding set 79 of fragility functions. Damage descriptions are commonly reduced to ordered, collectively 80 exhaustive and mutually exclusive damage states (e.g., no, slight, moderate, heavy and com- 81 plete damage), which are numerically encoded as ds∈ {0,1, ..., c}. To ease the notation 82 in the following definitions, we consider a single building class and assume that all damage 83 observations for buildings of that class are treated equally. 84 The log-normal cumulative distribution is the most common functional form for fragility 85 functions 86 P(DS≥ds|im) = Φln(im/θ ds) βds = Φlnim−lnθds βds , (1) where Φ(·)is the standard normal cumulative distribution function, θdsdenotes the median 87 IMwhich causes the structure to reach or exceed ds, and βdsis the standard deviation of 88 thelnIM, here referred to as the dispersion parameter. From the fragility functions, one 89 obtains the probability mass function, P(DS=ds|im) =p(ds|im), as 90 p(ds|im) =  1−P(DS≥ds+ 1|im), ifds= 0 P(DS≥ds|im), ifds=c P(DS≥ds|im)−P(DS≥ds+ 1|im),otherwise .(2) 4 To prevent negative probabilities in Eq. (2), we assume an identical dispersion parameter 91 βfor all damage states. This avoids a crossing of the fragility functions and is a standard 92 assumption taken in analytical and empirical fragility modeling for ordered damage states 93 (e.g., Martins and Silva, 2020; Rosti et al., 2021). For c+ 1damage states, we thus have c 94 increasing parameters θ1< θ 2<···< θcand one common dispersion parameter β. 95 A second fragility function definition, equivalent to Eq. (1), is obtained from the cu- 96 mulative probit model (Lallemant et al., 2015). By introducing the threshold parameters 97 ηds=β−1lnθds, this alternative definition reads as 98 P(DS≥ds|im) = Φ β−1lnim−ηds . (3) FollowingNguyenandLallemant(2022), weusethere-parametrizeddefinitionofEq.(3)with 99 the transformed parameters ηdsstill being of increasing order. To account for this ordering 100 in the estimation algorithms, we introduce positive parameters δds=ηds−ηds−1for states 101 ds > 1. The parameters of interest are collectively denoted as vector ϑ= [β, η 1, δ2, . . . , δ c]⊤. 102 Finally, we denote the fragility model in terms of the probability mass function p(ds|im,ϑ), 103 where we explicitly condition on parameter vector ϑ. 104 The use of log-normal fragility functions – leading to the cumulative probit model – dates 105 back to earthquake risk assessments of critical infrastructure components (e.g., Kennedy 106 et al., 1980) and is well-established in disaster risk analysis (e.g., Kircher et al., 1997) as 107 well as in performance-based earthquake engineering (e.g., Porter et al., 2007). Yet, the 108 herein presented estimation algorithms can also be used for alternative models, such as the 109 cumulative logistic model employed by Basöz and Kiremidjian (1998). 110 5 ESTIMATING FRAGILITY FUNCTION PARAMETERS 111 Figure 1 depicts the probabilistic model relating the variables involved in the estimation of ϑ 112 by considering a damage data set from nbuildings gathered after an earthquake with rupture 113 characteristics rup. The fragility model of Eq. (2) relates the observed damage states, 114 ds= [dsi, . . . , ds n]⊤, to the IMvalues at the locations of the surveyed buildings, imB= 115 [im1, . . . , im n]⊤. Knowing the latter would allow estimating ϑwith similar techniques used 116 in analytical fragility modeling (e.g., Baker, 2015). In most cases, however, we only know the 117 IMvalues at the locations of seismic stations, here denoted as imS= [imn+1, . . . , im n+m]⊤118 for the case of mstations, while their counterparts at the survey locations are uncertain. 119 … … Station sites Building damageRupture characteristics Parameters of fragility functionsIntensity measure s Building sites … Figure 1: Probabilistic model to estimate fragility function parameters using data from building damage surveys and seismic network stations. Observed and unobserved variables are indicated with gray and white circles, respectively. The arrows indicate directed rela- tionships between the variables. The horizontal lines illustrate that the intensity measures follow a joint multivariate distribution. To constrain the above-mentioned uncertainty, we evaluate the distribution of imBcondi- 120 tional on station data imSand rupture characteristics rupby following Engler et al. (2022), 121 6 who present the mathematical background for the most recent version of the United States 122 Geology Survey’s shake map system. Conditional on ruponly, the joint distribution of 123 imBandimSis assumed to be multivariate log-normal (Jayaram and Baker, 2008), with 124 parameters derived from empirical ground motion models (GMMs) and spatial correlation 125 models. Thus, the conditional distribution of interest also follows a multivariate log-normal 126 distribution, i.e., p(imB|imS,rup) =LN(µB|S,ΣBB|S). The distributional parameters, µB|S127 andΣBB|Sare computed via Equations 4 and 5 in Engler et al. (2022), also explained in 128 Appendix A of this study. 129 We next present two methodologies to estimate ϑ. The first, traditional, method assumes 130 fixed – or deterministic – values for imB. The second, proposed, method, treats both, ϑand 131 imB, as uncertain and follows a Bayesian approach for inference. 132 Maximum likelihood inference with fixed IM values 133 This estimation method replaces uncertain IMvalues at the survey sites with fixed, best- 134 estimate, values derived from the conditional distribution p(imB|imS,rup). For site i, the 135 best-estimate value corresponds to cimi= exp([ µB|S]i), which is the median IMconditional 136 on station data imS. As depicted in Figure 1, the joint distribution of building damage 137 conditional on the IMvalues at the building locations factorizes, i.e., p(ds|imB,ϑ) = 138 Qn i=1p(dsi|imi,ϑ). Assuming fixed IMvalues thus offers the advantage that the damage 139 observations of different buildings can be treated independently, which makes it an efficient 140 approach for large amounts of data. The parameters are then estimated by maximizing the 141 log-likelihood 142 ˆϑ= argmax ϑnX i=1lnp(dsi|cimi,ϑ), (4) such that the obtained parameter estimates maximize the probability of observing the sur- 143 veyed damage under the chosen best-estimate IMvalues. In the following, we refer to this 144 7 estimation method as the fixed IMapproach. 145 Figure 2 illustrates a one dimensional example with simulated IMobservations and dam- 146 age data from two stations and 500 buildings resulting from a magnitude 6 earthquake with 147 an epicenter at x= 0km. The assumed fragility parameter values are 0.65, 0.2g, and 0.65g 148 forβ,θ1andθ2, respectively, and the estimated values are 0.92, 0.25g, and 0.80g. From 149 Figure 2b, we observe the substantial uncertainty around the median of p(imB|imS,rup). 150 By considering fixed IMvalues, this approach does not account for deviations of the “true” 151 IMfrom this best-estimate. The estimated fragility functions, shown in Figure 2c, are flatter 152 than the “true” functions used to generate the data. Such behavior was also observed in the 153 simulated case studies by Ioannou et al. (2015) mentioned in the Introduction, and results 154 from an overestimation of the dispersion parameter β. This overestimation indicates that 155 IMuncertainty manifests as apparent uncertainty – or dispersion – in the fragility functions. 156 012ds damage survey dataa) 20 15 10 5 0 5 10 15 20 x [km]102 101 100im [g]b) median 90% CI true station 0.0 0.2 0.4 0.6 0.8 1.0 im [g]0.00.20.40.60.81.0P(DS≥ds|im)ds=1 ds=2c) fixed IM estimate true Figure 2: Fragility results from a one-dimensional example with simulated damage survey data. (a) Building damage states plotted versus their location. Damage states are simulated using the assumed “true” fragilities and IMvalues. (b) Simulated “true” IMvalues, with circles denoting the locations and IMvalues for two recording stations. The dashed line in- dicates the predicted median IMvalues and the shaded interval indicates the 90% credibility interval (CI) of the distribution conditional on data from the two stations. (c) The “true” fragility functions and the ones estimated with the fixed IMapproach. 8 Bayesian inference with uncertain IM values 157 To account for ground motion uncertainty, we take a Bayesian approach treating both the 158 fragility function parameters, ϑ, and the IMvalues at the survey sites, imB, as realizations 159 of random variables. The model shown in Figure 1 serves as the prior generative model 160 which gives rise to the observed data. Based on this model, we express the joint posterior 161 distribution of imBandϑconditional on the triplet of damage survey data, ds, station data, 162 imS, and rupture characteristics, rup, as 163 p(ϑ,imB|ds,imS,rup) = nY i=1p(dsi|imi,ϑ)! p(imB|imS,rup)p(ϑ) p(ds|imS,rup), (5) where the denominator, p(ds|imS,rup), is the marginal likelihood, i.e., the probability that 164 the prior model assigns to the observed damage data conditional on station data and rupture 165 characteristics, and p(ϑ)is the prior distribution of the fragility function parameters. For 166 each parameter we assign weakly informative priors, stated in Table 1, and assume that the 167 parameters are a-priori independent. These distributions are chosen such that the result- 168 ing fragility functions provide sufficient adaptability to cover building classes with varying 169 susceptibility to earthquake damage, as illustrated in Appendix B. 170 Table 1: Prior distributions for the parameters ϑand corresponding means and the 5% and 95% quantiles ( q0.05andq0.95). Parameter Domain Distribution Mean q0.05q0.95 β R+p(β) = InvGamma( α= 2.5, β= 1.0)0.77 0.20 2.00 η1 R p(η1) =N(µ=−1.61, σ= 1.5) -1.61 -4.08 0.86 δds∀ds > 1R+p(δds) = Gamma( α= 1.5, β= 2.5) 0.6 0.07 1.56 The posterior distribution in Eq. (5) is not analytically tractable, thus we use Markov 171 Chain Monte Carlo (MCMC) sampling for Bayesian inference. MCMC is a computational 172 techniqueusedtoapproximatecomplexprobabilitydistributions. Itcombinesprinciplesfrom 173 9 MarkovchainsandMonteCarlomethodstoefficientlysamplefromatargetdistributionwhen 174 direct sampling is not feasible. This makes it particularly valuable for Bayesian statistics, 175 where we aim to draw samples from a certain target posterior. Here, we employ the No-U- 176 Turn Sampler (Hoffman and Gelman, 2014) as implemented in the software library NumPyro 177 (Phan et al., 2019). We obtain 3,000 sample sets (ϑ,imB)from the posterior through four 178 independent MCMC chains with 1,000 warm-up steps each. The convergence of the chains is 179 tested by following the recommendations of Vehtari et al. (2021) summarized in Appendix C. 180 The implementation details can be found in the provided code repository1. 181 Figure 3 illustrates the posterior samples of (ϑ,imB), obtained from MCMC, for the 182 same one-dimensional example used in Figure 2. From Figure 3b, we observe how the 183 consideration of the damage survey data allows for an improved estimate on imBcompared 184 tothecaseinFigure2bthatusessolelytheavailablestationdata. Figure3cplotsthefragility 185 functions obtained from the mean of the sampled parameters (dashed lines, corresponding 186 to 0.65, 0.19g, and 0.61g for β,θ1andθ2, respectively), which closely match the “true” data- 187 generating functions. The shaded areas indicate the 90% credibility interval (CI). The latter 188 is the difference between the 95% and 5% quantiles obtained at each imlevel for all functions 189 computed from the posterior parameter samples. The limited amount of data considered in 190 this simplified example results in relatively large CIs, especially for ds= 2. 191 Accounting for different building characteristics 192 So far, we have assumed identical fragility functions for all considered buildings. To reflect 193 the heterogeneity in realistic building stocks, the model should account for the differences 194 in structural types, geometry, and other building characteristics. For this purpose, most 195 current regional earthquake risk models use multiple fragility models, each defined for a 196 specific combination of building characteristics. These combinations are here referred to as 197 1www.doi.org/10.5281/zenodo.10074233 10 012ds damage survey dataa) 20 15 10 5 0 5 10 15 20 x [km]102 101 100im [g]b) median 90% CI true station 0.0 0.2 0.4 0.6 0.8 1.0 im [g]0.00.20.40.60.81.0P(DS≥ds|im)ds=1 ds=2c) Bayesian estimate trueFigure 3: Fragility results from the Bayesian approach for the example from Figure 2: (a) Simulated building damage states plotted versus their location. (b) Simulated “true” IM values, with circles denoting the locations and IMvalues for two recording stations. The dashed line and the shaded area indicate the median and the 90% CI of the posterior IM distribution conditional on station and damage data. (c) The “true” and the estimated fragility functions. The dashed lines show the functions obtained from the mean posterior parameters, and the shaded area indicates the 90% CI. building classes and denoted with variable bc∈ BC, where BCis the set of all considered 198 building classes. 199 The fragility functions for different bcare commonly assumed to be of the same functional 200 form, i.e., the model specified in Eq. (3), but with differing parameters. Thus, the total set of 201 parametersis ϑ={ϑbc|bc∈ BC}. Whenusingthefixed IMapproach, theparameters ϑbcare 202 estimated independently for each bc. Thus, it may be necessary to use fewer building classes 203 such that a sufficiently large number of observations per class is available. The proposed 204 Bayesian approach jointly estimates the posterior of IMvalues and all parameters ϑ. This 205 results in a more complex estimation step, but has the advantage that even classes with 206 few observations profit from the reduced uncertainty in the IMposterior obtained from all 207 observations. 208 Information on building characteristics is gathered during the building inspection for 209 11 post-earthquake damage assessment. The herein presented case study makes use of existing 210 studies that defined a set of building classes, BC. In the general case, however, the choice 211 ofBCand the fragility modeling are inherently connected. It often requires an iterative 212 procedure to check which BCperforms best in distinguishing susceptibility of a building to 213 earthquake damage, while still providing enough observations per building class. 214 CASE STUDY: THE 2009 L’AQUILA (ITALY) EARTHQUAKE 215 Thiscasestudyestimatesfragilityfunctionparametersusingdatafromthe2009M6L’Aquila 216 earthquake. The objectives are twofold: first, to apply the proposed Bayesian estimation 217 method on a real data set and second, to examine the role of IMuncertainty by comparing 218 the Bayesian estimation results to results from the traditional, fixed IMapproach. The 219 data set contains results of visual post-earthquake inspections of residential buildings and 220 was published by the Italian civil protection department (Dolce et al., 2019). As such, it 221 providesinformationonthesustaineddamageanddescriptionsofthebuildingcharacteristics 222 (e.g., the geometry and the structural type). This survey data is combined with rupture 223 characteristics and station data from the engineering strong motion database (Luzi et al., 224 2020). Figure 4a shows the surface projection of the rupture and the stations located within 225 the region of interest together with soil conditions in terms of the 30 m time-averaged shear- 226 wave velocity, vs30, obtained from Mori et al. (2020). 227 The considered data set consists of 56,400 masonry buildings distributed over 62 munic- 228 ipalities, as shown in Figure 4b. The damage states are derived from the raw inspection 229 results based on Table 1 in Rosti et al. (2021) and consist of six categories: no ( ds= 0), 230 negligible to slight ( ds= 1), moderate ( ds= 2), severe ( ds= 3), very heavy damage ( ds= 4) 231 and collapse ( ds= 5). In accordance with the national Italian risk assessment (Dolce et al., 232 2021), we categorize the masonry buildings into three structural types, A, B and C1, with 233 12 Figure 4: Overview of the L’Aquila case study region: a) seismic network stations, projected rupture surface, and soil conditions as measured via vs,30; b) number of buildings per mu- nicipality in the considered data set. decreasingsusceptibilitytoearthquake-induceddamage. Thelinkbetweenspecificattributes 234 in the damage survey and the three structural types are described in Table 10 in Dolce et al. 235 (2019). Finally, fragility parameters are estimated for six building classes, derived from the 236 structural type and the number of storeys. Specifically, each type is further divided into 237 low-rise structures (less than two storeys) and mid-/high-rise structures (more than three 238 storeys), denoted as L and MH, respectively. We estimate six parameters for each building 239 class, leading to a total of 36 parameters. 240 The above-described data set was obtained after several pre-processing steps of the raw 241 post-earthquake survey results. One often-encountered challenge in empirical fragility mod- 242 elling is the under-representation of undamaged buildings. This occurs because inspections 243 may focus preferentially on damaged buildings. To alleviate this issue, we broadly follow the 244 13 approach taken by Rosti et al. (2021). For each municipality, the number of buildings in the 245 database is compared to statistics from the 2001 national census (ISTAT, 2001). Municipal- 246 ities where the available data set covers at least 90% of the number of buildings indicated 247 in the census statistics are considered as completely surveyed. For municipalities where this 248 ratio is less than 10%, it is assumed that only damaged buildings were inspected, and the 249 remaining buildings are added to the inspection data set assuming that they are undamaged. 250 Inspection results from municipalities with ratios between 10% and 90% are not considered, 251 due to the challenge of inferring the condition of non-inspected buildings in those cases. 252 Missing fields in the original inspection data set and the added undamaged buildings 253 result in some buildings having unknown spatial coordinates, number of storeys, or struc- 254 tural type. Buildings with identical or missing coordinates (e.g., the added buildings) are 255 re-distributed based on the European settlement map (Corbane and Sabo, 2019). Missing 256 entries on storey number and structural type are imputed based on census and survey data. 257 For this illustrative case study, we consider the derived inputs as fixed, but we could also 258 propagate such input uncertainty within the proposed Bayesian framework. Of the 56,400 259 masonry buildings in the final data set, 35,700 (63%) stem from the post-earthquake inspec- 260 tion data set and 20,700 (27%) undamaged buildings are added based on census data. For 261 each considered building class, Table 2 shows the number of buildings per damage state. 262 Table 2: Number of buildings of a certain building class and damage state for the compiled damage survey data set. building class ds= 0 ds= 1 ds= 2 ds= 3 ds= 4 ds= 5 A-L 8,915 2,771 1,219 1,855 2,059 1,570 A-MH 4,633 1,895 810 1,138 1,453 874 B-L 8,763 1,725 494 570 491 352 B-MH 4,871 1,263 377 430 444 290 C1-L 3,425 542 111 107 115 60 C1-MH 2,077 400 93 97 66 55 14 The pre-processed damage data set, summarized in Table 2, is then used to estimate 263 fragility functions using two approaches: the fixed IMapproach and the proposed Bayesian, 264 MCMC-based,approach. Thelatterusestheentireconditionaldistribution p(imB|imS,rup), 265 which requires the storage of the (n×n)-dimensional matrix ΣBB|S. Due to this high memory 266 demand, the size of the considered data set has to be limited. The inference is run on an 267 NVIDIA V100 tensor core GPU with 25 GB of RAM that allows the consideration of data 268 from n=14,000 buildings, corresponding to one fourth of the data set. The data subset is 269 developed by randomly sampling 25% of the buildings from each building class/damage state 270 combination shown in Table 2. For the computationally less expensive fixed IMapproach, 271 we use the entire data set. 272 We first estimate fragility functions using SA(0.3s)as the IM. Thus, the vectors imBand 273 imSare spectral accelerations at sites of surveyed buildings and network stations, respec- 274 tively. To compute the conditional distribution p(imB|imS,rup), we use the OpenQuake 275 (Silvaetal.,2014)implementationoftheBindietal.(2011)GMMandthe spatialcorrelation 276 model of Esposito and Iervolino (2012). Figure 5a plots the median spectral accelerations 277 conditional on station data imS. Recall that this median value is used to estimate param- 278 etersϑwith the fixed IMapproach. The Bayesian approach, on the other hand, provides 279 samples from the joint posterior of imBandϑ. The median of the posterior imBsamples is 280 plotted in Figure 5b, where the locations of the considered buildings are also indicated. We 281 observe that the resulting spatial IMpattern is more refined compared to the original map 282 in panel a. 283 Figure 6 plots the fragility functions of two building classes, A-MH and C1-L, as obtained 284 with the fixed IMapproach and the proposed Bayesian approach. Similar to the one- 285 dimensional example in Figure 2, we observe that the fixed IMbased parameter estimates 286 lead to flatter fragility functions compared to those obtained from the Bayesian approach, 287 providing higher damage probabilities for low IMvalues and lower probabilities for large 288 15 Figure 5: The spatial distribution of median SA(0.3s)values conditional on seismic net- work recordings (a) and additionally conditioned on damage survey data via the proposed Bayesian approach (b). Network stations and locations of surveyed buildings are indicated with triangles and points, respectively. IMvalues. For example, considering building class A-MH in Figure 6a, and a spectral 289 acceleration SA(0.3s) = 0 .25g, the fragility functions obtained from the fixed IMapproach 290 predict a 30% probability of reaching or exceeding ds= 4, while this probability decreases to 291 5% if the mean of the posterior samples is used. For the Bayesian approach, we also note the 292 decrease in the 90% CI compared to the one-dimensional example in Figure 3. This is due 293 to the larger size of the considered damage survey data set, combined with a larger number 294 of seismic stations. 295 Table 3 provides the estimated fragility function parameters for the six considered build- 296 ing classes. The fixed IMapproach leads to substantially larger estimates of the dispersion 297 parameter β, explaining the flatter fragility functions in Figure 6. By fixing IMvalues to 298 the best-estimate values plotted in Figure 5a, the variability of the actual, unobserved, IM 299 values around these best-estimates manifests in apparent dispersion in the estimated fragility 300 functions, as reflected by the large βvalues. Both estimation approaches capture the in- 301 16 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]0.00.20.40.60.81.0P(DS≥ds|sa)A-MHa) fixed IM Bayesian 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]C1-Lb)Figure 6: Fragility functions for building class A-MH (a) and C1-L (b) estimated from the damage survey data using the fixed IMand the Bayesian approaches. The five functions correspond to ds∈ {1,2, . . . , 5}from the top left to the bottom right. Functions from the Bayesian approach are illustrated using the mean posterior samples (solid lines) and the 90% CI (shaded area). creasing damage susceptibility expected from the employed building categorization scheme, 302 i.e., buildings of type A having lower θdsvalues and thus being more fragile than buildings 303 of type C1. Furthermore, low-rise buildings (e.g., A-L) tend to be less fragile than their 304 counterparts with more than two storeys (e.g., A-MH). 305 Table 3: Fragility function parameters for the six considered building classes estimated with the fixed IMapproach (FA) and the proposed Bayesian approach (BA). The stated Bayesian parameter estimates correspond to the mean of the posterior samples. Class β[-] θ1[g] θ2[g] θ3[g] θ4[g] θ5[g] FA BA FA BA FA BA FA BA FA BA FA BA A-L 1.09 0.68 0.18 0.30 0.31 0.46 0.39 0.55 0.60 0.73 1.13 1.13 A-MH 0.96 0.59 0.14 0.25 0.25 0.41 0.32 0.50 0.48 0.66 0.96 1.05 B-L 1.20 0.81 0.39 0.47 0.78 0.82 1.01 0.99 1.50 1.33 2.65 2.00 B-MH 1.20 0.75 0.29 0.37 0.60 0.65 0.78 0.79 1.15 1.05 2.19 1.62 C1-L 1.35 0.82 0.62 0.63 1.47 1.14 1.93 1.37 2.76 1.74 5.43 2.73 C1-MH 1.17 0.75 0.46 0.52 1.03 0.93 1.35 1.13 2.01 1.49 3.16 2.01 17 Results using alternate IMs 306 Only few empirical fragility functions are derived for spectral accelerations, while the ma- 307 jority are based on PGA. The use of SA(0.3s)above was motivated by the fact that it is 308 seen as a more efficient predictor of building response. Given this assumption, we would 309 expect smaller dispersion values, β, when using spectral accelerations. Figure 7 shows the β 310 values estimated for both IMs using the fixed IMand the Bayesian approaches. Indeed, the 311 βvalues for SA(0.3s)are slightly smaller for both approaches. The small difference may be 312 explained by the fact that we focus on masonry buildings that are almost all less than four 313 storeys tall, thus representing relatively heavy and non-ductile structures for which PGA 314 can be an adequate fragility predictor. For the remainder of this study, we continue with 315 SA(0.3s)as the IMof interest. 316 A-L A-MH B-L B-MH C1-L C1-MH0.00.20.40.60.81.01.21.41.6dispersion βfixed PGAfixed SA(0.3s)Bayesian PGA Bayesian SA(0.3s) Figure 7: Comparison of estimated βvalues from the fixed IMand the Bayesian approach forPGAandSA(0.3s). 18 Comparison for different GMMs 317 Sofar, weusedtheBindietal.(2011)GMMandtheEspositoandIervolino(2012)correlation 318 modeltoestablishthedistributionof IMvaluesconditionalonstationdata(illustratedbyits 319 median in Figure 5a). These models are calibrated on Italian ground motion data and used 320 in the official Italian shake map system (Michelini et al., 2020). We denote this combination 321 as GMM1. Here, we will compare the parameter estimates obtained with GMM1 to the 322 ones resulting from a second combination, named GMM2, that uses the Chiou and Youngs 323 (2014) GMM and the correlation model of Bodenmann et al. (2023). Both GMM2 models 324 are estimated on the same global ground motion data set. The change of input motion 325 models affects only the imBdistribution in Figure 1, and we use the alternate distribution 326 to then re-estimate the fragility parameters, ϑ. 327 Figure 8 shows the fragility functions for A-L buildings derived using GMM1 and GMM2, 328 respectively. For better readability, only the mean Bayesian estimates are shown. From 329 Figure 8b we observe that the functions obtained from the Bayesian approach are similar for 330 both GMMs, while Figure 8a indicates a substantial sensitivity on the underlying GMM for 331 the fixed IMapproach. In the latter case, the parameter estimates based on GMM2 indicate 332 less fragile buildings. The reasons for these differences are explained below. 333 The alternate models within GMM1 and GMM2 predict differing SA(0.3s)values for the 334 rupture due to differences in functional forms, dependent variables, and the data used in 335 fitting the models’ parameters. To illustrate these differences, we compare the SA(0.3s)pre- 336 dictions along the horizontal transect illustrated in Figure 9a. Figure 9b plots the predictive 337 distributions in terms of the median and the 90% CI, prior to taking into account station 338 recordings and damage survey data. For sites located within the rupture surface projection, 339 GMM2 predicts larger SA(0.3s)values than GMM1. These differences also persist after 340 taking into account station recordings, as visible from the corresponding distributions in 341 19 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]0.00.20.40.60.81.0P(DS≥ds|sa)A-Lfixed IM a) GMM1 GMM2 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]A-LBayesian b) GMM1 GMM1Figure 8: Comparison of fragility functions for building class A-L that are estimated using two different combinations of ground motion and spatial correlation models (GMM1 and GMM2) through the fixed IMapproach (a) and the proposed Bayesian approach (b). Figure 9c. Recall that the fixed IMapproach uses the median values indicated in Figure 9c 342 to estimate the fragility functions. This explains the differences in the estimated functions 343 discussed above and shown in Figure 8a. 344 Finally, Figure 9d plots the distributions as obtained after processing the damage data 345 through the proposed Bayesian approach. For locations in the vicinity of surveyed buildings, 346 we observe similar distributions from both GMMs. By performing joint Bayesian inference 347 over SA(0.3s)values and fragility parameters, the posterior SA(0.3s)distributions tend to 348 converge regardless of whether GMM1 or GMM2 is used to establish the prior distribution 349 ofSA(0.3s)values. This explains the similarity in the estimated fragility functions observed 350 in Figure 8b, and indicates that the proposed Bayesian approach is more robust than the 351 traditional, fixed IM, approach with respect to such modelling assumptions. 352 Effect of data sub-sampling 353 The proposed Bayesian approach is computationally expensive, and thus only one fourth 354 of the entire data set is used for parameter estimation. This is a drawback compared to 355 20 42°20'N42°25'N station transect survey10 kma) 0.020.050.10.20.51.02.0SA(0.3s) [g] 0.56g0.99gb)GMM1 GMM2 0.020.050.10.20.51.02.0SA(0.3s) [g] 0.44g0.93gc) 13°15'E 13°30'E 13°45'E 14°00'E0.020.050.10.20.51.02.0SA(0.3s) [g] 0.59g0.67gd)Figure 9: Comparison of SA(0.3s)median values and 90% CIs obtained from GMM1 and GMM2, andplottedalongahorizontaltransectthroughtheregion: (a)Stationsandsurveyed buildings in the vicinity of the chosen transect. (b) GMM predictions conditional on rupture characteristics; (c) after additional conditioning on station data; and (d) after additional conditioning on damage data using the Bayesian approach. The labeled points in (b to d) are example values for a numerical comparison. the fixed IMapproach that handles the entire data set with no computational difficulties, 356 and raises the question of how the Bayesian parameter estimates change for different sub- 357 sets. Figure 10 shows βvalues estimated with four randomly chosen sub-samples of the 358 entire data set. The four estimated values for each building class are within the posterior 359 distribution CIs, and with differences that are much less than the difference relative to the 360 fixed IMvalues, indicating that the proposed approach produces stable outputs even with 361 random sub-sampling. This is perhaps not surprising given the large number of observations 362 that remain within a given sub-sample. Recall that we randomly selected 25% of the data 363 21 points from each building class/damage state combination shown in Table 2. The use of 364 more sophisticated sub-sampling techniques is likely to further reduce the sensitivity of the 365 estimated parameters. 366 A-L A-MH B-L B-MH C1-L C1-MH0.00.20.40.60.81.01.21.41.6dispersion βfixed IM Bayesian (4 random subsamples) Figure 10: Comparison of estimated βvalues from the fixed IMapproach and for MCMC estimated from four random subsamples of the entire data set. In a separate analysis, we compare the performance of both estimation methods in recov- 367 ering a set of known, “true”, fragility functions while using the same sub-sampling technique 368 as above. For this purpose, we replace the observed damage states in the L’Aquila data set 369 by simulated values. These are obtained by first drawing one sample of imBfrom the estab- 370 lished conditional distribution p(imB|imS,rup)and then using a pre-defined set of fragility 371 functions to sample a damage state for each building. Figure 11 compares the estimated 372 fragility functions of two building classes, A-L and B-MH, to the “true” data-generating 373 functions. For better readability, we only plot the functions for ds∈ {1,2,4}. The functions 374 obtained from the fixed IMapproach are flatter than the “true” functions, while the Bayesian 375 approach leads to functions that are in good agreement, despite only one fourth of the data 376 set being used. 377 22 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]0.00.20.40.60.81.0P(DS≥ds|sa)A-La) fixed IM Bayesian true 0.00 0.25 0.50 0.75 1.00 1.25 1.50 SA(0.3s) [g]B-MHb)Figure 11: Comparison of fragility functions of building class A-L (a) and B-MH (b) es- timated from a simulated version of the L’Aquila damage survey data set to the “true” data-generating functions. For better readability, only functions for ds∈ {1,2,4}are plot- ted (from top left to bottom right). CONCLUSION 378 This study explored the role of IMuncertainty in the estimation of fragility function pa- 379 rameters from damage survey data gathered after earthquake events, and proposed a novel 380 estimation method that can account for this uncertainty. Traditional estimation methods 381 do not allow to account for such uncertainty and instead assume fixed, deterministic, IM 382 values. In simulated case studies, traditional methods reportedly failed to recover the “true” 383 data-generating fragility functions. To address this shortcoming, we proposed a Bayesian 384 estimation approach that simultaneously considers IMvalues and fragility function parame- 385 ters as uncertain variables. We derived an expression for their joint distribution conditional 386 on the gathered damage data. For parameter estimation, we approximate this posterior 387 distribution through MCMC sampling. 388 We applied the proposed method in a case study using damage data gathered after the 389 2009L’Aquilaearthquakeevent. Thedatasetconsistsof56,400residentialmasonrybuildings 390 23 that are categorized into six building classes and for which overall damage is described with 391 six categories ranging from no damage to collapse. 392 We used both synthetic data and the L’Aquila data to compare the Bayesian parameter 393 estimates with estimates from a traditional, fixed IM, approach. The two methods produce 394 surprisingly different fragility functions, with those obtained from the fixed IMapproach 395 being flatter than their counterparts derived with the proposed Bayesian approach. The 396 flatter fragility functions predict higher probabilities of damage for low IMvalues and, con- 397 versely, lower probabilities for high IMvalues. Such behavior is observed across all damage 398 states and building classes. By fixing IMvalues to best-estimate values, the variability of 399 the actual, unobserved, IMvalues around these best-estimates manifests in an increased ap- 400 parent dispersion in the estimated fragility functions. The proposed Bayesian approach, on 401 the other hand, explicitly accounts for both variability sources and performs joint inference 402 onIMvalues and model parameters, which prevents such an inflation of fragility function 403 dispersion. 404 Additional tests demonstrated that the proposed Bayesian approach is less sensitive to 405 the choice of ground motion and spatial correlation models used to derive the distribution 406 ofIMconditional on station data (i.e., the shake map) than the traditionally used method, 407 which relies on fixed IMvalues. This reduced sensitivity is a significant advantage, as the IM 408 predictions from any given ground motion and spatial correlation model are not “correct”, 409 as must be assumed with the fixed IMmethod. The proposed method of jointly estimating 410 theIMvalues and fragility parameters allows for more robust estimates with respect to such 411 modeling assumptions. 412 TheproposedBayesianapproachcomeswithcomputationalchallenges, inparticularwith 413 respect to memory demand. The herein used implementation could perform inference for 414 14,000 data points. This number is sufficient for many damage data sets, but its application 415 to the L’Aquila data set required sub-sampling the data set. To assess the implications 416 24 of such sub-sampling we performed two analyses. First, we found that the parameters 417 estimated from multiple randomly chosen subsets were relatively stable across the different 418 subsets. Second, we replaced the observed damage with simulated values and found that 419 the Bayesian approach could recover the “true” data-generating fragility functions using only 420 a subset of the data. The fixed IMapproach, on the other hand, failed to do so even 421 when the entire data set was used for parameter estimation. While the proposed method 422 does require an additional analyst effort and more computational resources relative to the 423 traditional estimation method, the required effort is still vastly smaller than the time and 424 cost of inspections to collect the damage data, and so is a good investment given the extra 425 value it provides. 426 While the focus of the presented work lies on IMuncertainty, the proposed method can 427 be extended to account for uncertainties in other model inputs, such as missing database 428 entries. The method could also be adapted to fragility models that have different dispersion 429 parameters for each damage level and depend on building characteristics in a more complex 430 manner than the commonly used categorization into pre-defined building classes. Despite 431 these opportunities for further research, the novel estimation method can already be of 432 substantial benefit for analysts interested in improved empirical fragility modeling. For 433 the above reasons, we strongly believe that the method’s presented comparative advantages 434 justify its increased computational complexity. 435 Software tool To facilitate the application of the proposed Bayesian parameter estimation 436 approach, we developed a user-friendly Python-based software tool, which is available at 437 www.doi.org/10.5281/zenodo.10074233. The tool comes with several tutorials and example 438 data sets that further explain the features and implementations of the proposed estimation 439 procedure. Importantly, it contains the code to reproduce the herein presented results and 440 figures. 441 25 Acknowledgements The first author gratefully acknowledges support from the ETH Risk 442 Center (grant 2018-FE-213). 443 APPENDIX 444 A DISTRIBUTION OF IM VALUES CONDITIONAL ON STATION DATA 445 This appendix provides further details on the approach of Engler et al. (2022) that was used 446 in this study to compute the distribution of IMvalues at locations of surveyed buildings, 447 imB, conditional on IMvalues at locations of seismic network stations, imS, and rupture 448 characteristics rup. For this purpose, we first look at the distribution of IMvalues at a 449 generic set of spatially distributed sites conditional on rupand in absence of station data. 450 These values are assumed to follow a multivariate log-normal distribution (Jayaram and 451 Baker, 2008), in other words the lnIMvalues follow a multivariate normal distribution with 452 mean vector µand covariance matrix Σ. For each site i, empirical GMMs allow computing 453 the mean, µi, oflnIM iand the standard deviations, τiandϕi, of the between- and within- 454 event residual, respectively. The between-event residual is constant for any two sites iand 455 j, while the within-event residuals are partially correlated with coefficient ρij. The latter 456 is computed with a spatial correlation model. The entries of the covariance matrix are 457 [Σ]ij=τiτj+ϕiϕjρij, while the entries of the mean vector µfollow from the predicted mean 458 values µiobtained from the GMM. Thus, the joint distribution of IMvalues at survey and 459 station sites is a multivariate log-normal distribution expressed as 460 p  imB imS |rup =LN  µB µS ; ΣBBΣBS ΣSBΣSS  , (A1) where µBandΣBBdenote the mean vectors and covariance matrices of lnIMvalues at the locations of the surveyed buildings, µSandΣSStheir counterparts at the station locations, 26 andΣBSthe covariance matrix between the lnIMvalues at the buildings and the stations. Then, the parameters of p(imB|imS,rup)follow as µB|S=µB+ΣBSΣ−1 SS(lnimS−µS) (A2a) ΣBB|S=ΣBB−ΣBSΣ−1 SSΣSB, (A2b) where µB|SandΣBB|Sdenote the mean vector and covariance matrix of lnIMvalues at the 461 surveyed buildings conditional on station data imS. 462 B PRIOR DISTRIBUTION OF FRAGILITY FUNCTION PARAMETERS 463 Figure B1 illustrates the prior distributions for the fragility function parameters ϑstated 464 in Table 1 and plots the fragility functions obtained from these priors using an example of 465 three damage states, where ϑ= (β, η 1, δ2). For each ϑ, the fragility functions are obtained 466 by computing η2=η1+δ2, and Eq. (3). Alternatively, we may compute the median IM 467 which causes the structure to reach or exceed dsasθds= exp( βηds), and then use Eq. (1). 468 C MCMC CONVERGENCE DIAGNOSTICS 469 This appendix provides a brief summary on the conducted tests to verify the convergence of 470 the target distribution samples obtained from Markov Chain Monte Carlo (MCMC) simula- 471 tions. Sequences of iterative simulations can fail to converge, either because the chains are in 472 different parts of the target distribution or because the chains have not attained stationarity. 473 Thus, it is crucial to run multiple independent chains for convergence diagnostics. As a first 474 visual convergence assessment, one may compare the distribution of posterior samples from 475 each chain, as well as their sample traces. In the presented case study, we performed MCMC 476 over 14,036 variables (36 damage model parameters and the IMs at 14,000 sites). Performing 477 a visual assessment for each single variable is not practicable for such high-dimensional exer- 478 27 0 1 2 β0.00.51.01.5probability densitya) 4 2 0 η1b) 0 1 2 δ2c) 0.0 0.5 1.0 1.5 im [g]0.00.20.40.60.81.0P(DS≥ds|im)ds=1 d) 0.0 0.5 1.0 1.5 im [g]ds=2 e) sample 1 sample 2 sample 3 meanFigure B1: Illustration of the prior distributions for fragility function parameters β,η1,δ2 (a to c), and of the fragility functions obtained from the prior mean parameters, as well as three sampled parameter sets, for: d) ds= 1, and e) ds= 2. cises. In such cases, we analyze quantitative metrics that indicate poor convergence and flag 479 potential problems. The conducted diagnostics use information between and within chains 480 and are based on Vehtari et al. (2021) to which an interested reader is referred to for further 481 details. The diagnostic metrics are computed with the Phython-based software library ArviZ 482 (Kumar et al., 2019). 483 The arguably most popular metric is the potential scale reduction factor, ˆR, introduced 484 by Gelman and Rubin (1992), which compares the variance between multiple chains to the 485 variance within each chain. If the chains converge, the between-chain and within-chain vari- 486 ances should be identical, i.e., ˆR≈1. In addition to ˆR, we also compute the effective 487 sample size ( ESS) of a quantity of interest. In simple terms, ESSstates how many inde- 488 pendent draws contain the same amount of information as the dependent sample obtained 489 28 via MCMC. Both measures, ˆRandESS, are computed based on rank-normalized samples 490 as recommended by Vehtari et al. (2021). Based on the heuristics presented by the same 491 study, good convergence is achieved if ˆR≤1.01andESS>400. Consequently, we chose the 492 number of samples (750 per chain) and warm-up steps (1,000 per chain) such that both of 493 these conditions are fulfilled. 494