Thisstudyintroducesacomputationalframeworkforeﬃcientandaccurateseismicfragilityanalysisbased on a combination of artiﬁcial ground motion modeling, polynomial-chaos-based global sensitivity analysis, andhierarchicalkrigingsurrogatemodeling. TheframeworkfollowsthephilosophyofthePerformance-Based EarthquakeEngineeringPEERapproach, wherethefragilityanalysisisdecoupledfromhazardanalysis. This study addresses three criticalities that are present in the current practice. Namely, reduced size of hazard- consistent size-speciﬁc ensembles of seismic records, validation of structural simulators against large-scale experiments, high computational cost for accurate fragility estimates. The eﬀectiveness of the proposed framework is demonstrated for the Rio Torto Bridge, recently tested using hybrid simulation within the RETRO project. Keywords: Fragility analysis ; surrogate modeling ; hierarchical kriging ; polynomial chaos expansion ; global sensitivity analysis ; artiﬁcial ground motion. 1 Introduction The PEER Performance-Based Earthquake Engineering (PBEE) framework is widely used to perform seismic risk analysis. The framework is a simple statement of the total probability theorem, which combines the output of Probabilistic Seismic Hazard Analysis (PSHA) with fragility, damage, and loss analysis. A critical part of the framework is the coupling between hazard and structural modeling, namely, fragility analysis. It comes as no surprise that it developed into a stand-alone research ﬁeld. Given a structural system of interest, the most recent advancement in fragility analysis requires the selection of a set of suitable ground motions to be used as input for non-linear structural analysis. Then, the statistics of a Quantity of Interest (QoI), which in general is denoted as Engineering Demand Parameter (EDP), are derived as a function of a given Intensity Measure (IM), or a vector of IMs. Speciﬁcally, these conditional EDP-IM relations are named fragility functions. Following this line, in recent years, a lot of research has been devoted to creating computational fragility frameworks, including the selection and the use of real seismic records and non-linear time history analysis. An incomplete list of studies in this topic includes [1, 2, 3, 4]. However, within this perimeter, a number of open issues have been raised by the opinion paper of Bradley [5] and further elaborated in the state-of-the-art review of Silva and co-authors [6]. Departing from these studies, this paper aims to address three criticalities: 1. In recent years, there have been signiﬁcant advancements in ground motion selection procedures, which have partially solved the problem of deﬁning a suitable set of seismic records for structural analysis. However, although these selections are suﬃcient to describe important characteristics of the hazard for a given structure of interest, the number of available records is usually insuﬃcient to compute accurate estimates of the QoIs statistics. ∗Corresponding Author: Dr. Giuseppe Abbiati, abbiati@eng.au.dk , ORCID: 0000-0002-5048-8505 . 1 DRAFT2. The validation of structural simulators and the quantiﬁcation of the associated model uncertainty against representative experiments—despite being critical—is often overlooked. This is due to practical reasons. Structural testing is usually not aﬀordable beyond the component level due to the limited capacity of experimental facilities. Also, the budget allocated for structural testing is typically suﬃcient for running a small number of destructive experimental tests, which cannot cover the entire range of IMs for the fragility analysis. 3. The computation of fragility analysis requires a large number of time history analyses. It follows that computational cost of structural simulators limit de facto the total number of simulations. This is a classical problem in simulation-based Uncertainty Quantiﬁcation (UQ), which is tackled by replacing the computationally expensive simulator with an equivalent (with respect to, w.r.t., the QoI) response surface a.k.a. surrogate model. However, deﬁning surrogate models for non-linear dynamic analysis is challenging and still not completely addressed by the current state of the art of fragility computation. As a consequence, the accuracy of the fragility is limited by available computational resources. This study introduces a computational framework for eﬃcient and accurate seismic fragility analysis based on a combination of artiﬁcial ground motion modeling and surrogate modeling. Moreover, the goal is to integrate the know-how developed in ground motion selection and to follow the original spirit of the PEER-PBEE framework to decouple the task of hazard and fragility analysis. Speciﬁcally, given a suitable set of real seismic time series, we deﬁne an artiﬁcial ground motion model equipped with engineering-meaningful parameters. The selected real seismic time series are thought of as realizations of the ground motion model, using the same philosophy introduced by [7, 8]. Unlike their work, this framework uses the real seismic time series to calibrate the parameters of the artiﬁcial ground motion model and the associated uncertainties. The deﬁnition of such a hazard model poses the basis for using a surrogate-based UQ forward analysis. The building blocks of the computational framework are: the artiﬁcial ground motion model, an expensive- to-evaluate High-Fidelity (HF) structural simulator of the reference structure validated against experiments (e.g., Hybrid Simulation (HS) in the proposed application [9]), and a corresponding cheaper Low-Fidelity (LF) structural simulator (e.g., obtained via dynamic substructuring of the HF simulator in this study [10]). The seismic vulnerability analysis follows these sequential steps: i) generation a family of artiﬁcial motions (compatible with the real selected ones); ii) computation of the Polynomial Chaos Expansion (PCE) of the LF simulator QoI prediction; iii) dimensionality reduction of the artiﬁcial ground motion parameters based on Global Sensitivity Analysis (GSA); iv) computation of a Hierarchical Kriging (HK) surrogate that fuses LF and HF predictions of the QoI; v) the fragility analysis is derived as a natural byproduct of the framework. Fragility analysis and, more broadly, seismic risk assessment based on stochastic simulations is not new, and an incomplete list of studies includes [11, 12, 13, 14]. Moreover, to the best of our knowledge, the ﬁrst study using also kriging surrogate modeling is [15] and more recently [16]. Diﬀerent form the previous contributions, this paper draws a direct link between ground motion selection and the calibration of the stochastic hazard model; moreover, it introduces HK to combine predictions of simulators with diﬀerent degrees of ﬁdelity. In essence, the artiﬁcial ground motion model is used to build a large ensemble of time series, which would not be possible to compose using records of real seismic events. Therefore, this study aims to bridge classical seismic fragility analysis with the most recent UQ advancements, and not to solve the full seismic risk analysis (as done in [15] by fully employing the philosophy of [8]). In this context, this paper introduces a combination of PCE-based GSA and HK so that the computational burden associated with the fragility analysis is drastically reduced. The framework is applied to a real-world application consisting of a Reinforced Concrete (RC) bridge case study, which was recently tested within the RETRO project [17, 18]. LF and HF simulators are introduced, highlighting the crucial role of HS in enabling experimental validation of the bridge model [19]. Since the selected case study structure is insensitive to broadband ground motions, the artiﬁcial ground motion model incorporates a pulse-like component [20], and it is calibrated against a set hazard-consistent site-speciﬁc real records. The framework is applied to compute the fragility analysis related to the lateral drift of one of the twelve piers of the Rio Torto Bridge. This paper is organized as follows. Section 2 describes the proposed seismic fragility analysis framework. Section 3 introduces the Rio Torto Bridge case study. Section 4 describes both HF and LF structural simulators of the bridge. Section 5 describes the calibration of the artiﬁcial ground motion model using hazard-compatible site-speciﬁc seismic records. Section 6 summarizes the results of the seismic fragility analysis of the bridge. Section 7 includes limitations, future perspectives, and conclusions. For the sake of completeness, the entire machinery of surrogate modeling and, in particular, PCE-based GSA [21, 22] and HK [23] are reported as appendices providing a self-consistent notation. 2 DRAFT2 Seismic fragility analysis framework 2.1 Theoretical background LetA(t)denote the stochastic process representing the input excitation of the system. The goal is to represent the process in terms of an aleatory component (i.e., the inherent variability of the process) and an epistemic component (i.e., the parametric variability of the process). Accordingly, it is convenient to express A(t)as follows. A(t) =Ma(t;ZjXa); (1) whereMais a parametric stochastic model (or synthesis formula), Zis a set of random variables representing the aleatory variability, and Xaa is a set of model parameters. In this setting, the model parameters are considered as random variables. However, this uncertainty is epistemic (i.e., can be reduced). Let Ydenote the response QoI of a real system (e.g., the maximum displacement or internal force of a structural component), andMc(A(t)jXc)a suitable (w.r.t. the QoI) structural simulator, where Xcis a set of model parameters. It follows that Y, can be written as Y=Mc(Ma(t;ZjXa)jXc); (2) =Ms(ZjX); (3) whereMs:=McMaandX= [Xa;Xc]. Observe that while Mcis a deterministic structural simulator, Ms is a stochastic simulator. In fact, given a realization of the model parameters, x,Yis still a random variable with aleatory uncertainty determined by the propagation of Zjx. Using Eq.3, the Cumulative Distribution Function (CDF) of Y, denoted by FY(y), can be written as FY(y) =P(Ms(ZjX)y) =Z xZ zI(Ms(zjx)y)dFZ(z)dFX(x); (4) where I()denotes the indicator function, P()denotes probability, FZ(z),FX(x)denote the joint CDFs, and dFZ(z) =fZ(z)dz,dFX(x) =fX(x)dx, withfX(x)andfZ(z)being the joint probability densities. Observe thatZisstatisticallyindependentof X; moreover, fX(x) =fXc(xc)fXa(xa)underthe(reasonable)assumption that the parameters of the seismic input and the structural simulator are statistically independent. Next, the conditional CDF of Ygiven a realization xis simply FYjX(yjx) =P(Ms(ZjX=x)y) =Z zI(Ms(zjx)y)dFZ(z): (5) Given the distribution function FYjX, thep-quantile function returns the value ypsuch thatFYjX(ypjx) =p. AssumingFYjXis continuous and strictly monotonically increasing, yp=F1 YjX(pjx). However, since Xis a random variable, ypis also a random variable and can be written as Yp=F1 YjX(pjX) =Mp(X): (6) where, for convenience, we deﬁne Mp() :=F1 YjX(pj), andp2[0;1]. Therefore, Ypcan also be viewed as stochastic process indexed by p. Moreover,Mpis a deterministic model, since the aleatory variability has been marginalized in Eq.6. Given this general framework, a few observations need to be done for this speciﬁc study. First of all, A(t)is artiﬁcialgroundmotionmodel,and Xaisaphysicallymeaningfulsetofgroundmotionsparameters. Speciﬁcally, the model is a Gaussian ﬁltered white noise, therefore, Zis a set of standard normal random variables (more details in Section 5). Second, Mcis a structural simulator, and Xcis the set of model parameters. Finally, observe that despite its simplicity, Eq.6 involves a complex and highly non-trivial computational endeavor. In the following, this task is tackled using a non-intrusive UQ approach that assumes that Mc(and thereforeMs) is a black-box solver. 2.2 Computational solution In the context of earthquake engineering, the goal of this study is to determine the distribution of Ypfor diﬀerent values ofp. In general, the analytical solution of Eq.6 is not available; moreover, a classical Monte-Carlo-based 3 DRAFTsolution is computationally prohibitive. In fact, Mcis an "expensive" model and, therefore, also Mp. A natural alternative to drastically reduce the computational cost is to replace the original model with a surrogate model. Let ^Mpdenotes a surrogate model of Mp, then ^Yp=^Mp(X); (7) where ^YpYp. The statistics of ^Ypare then directly used for decision making. The surrogate model should be compatible with the artiﬁcial ground motion model, an HF validated structural simulator and the chosen QoI. Speciﬁcally, the proposed scheme builds on a combination of PCE-based GSA and HK, and it is here summarized in the following steps. 1.Model deﬁnition and validation. Deﬁne the QoI, a stochastic model of the input Ma(Section 5), and a set of structural simulators Mc(Section 4). In this study, c2fLF;HFgwherec=HFcorresponds to the expensive-to-evaluate HF simulator whereas c=LFto a computationally cheaper LF simulator (Section 2.2 provides further detail on this choices). The deﬁnition of the models reﬂects the following criteria: (a) The artiﬁcial ground motion model is deﬁned w.r.t. a local seismic catalog based on an existing suitable ground motion selection (e.g., [24, 25, 26]). In fact, this framework is designed to work complementary to any ground motion selection criteria, rather than to express a preference for speciﬁc selection criterion or to propose a new one. (b) The HF simulator successfully passed a comprehensive validation protocol. Validation entails a com- parison of model predictions with physical experiments that adequately reproduce the real structural behavior [27]. This study uses HF simulators validated using HS. This minimizes the cost of con- ducting realistic structural testing campaigns compared to shake table testing [9]. To be robust to generalization error, the HF simulator incorporates an accurate description of the physics of the real structure, which includes geometric and material properties. Therefore, it is ranked ﬁrst in terms of computational cost for a single evaluation amongst all structural simulators Mc. As a result, the HF simulator is expected to succeed in representing the structural response with predetermined accuracy within a given validation domain. (c) The LF simulator is a computationally cheaper proxy of its HF counterpart. In principle, there are no constraints on the types of allowed model simpliﬁcations. For example, one could decide using diﬀerent modeling paradigms for HF and LF simulators or 2D versus 3D geometrical representations. However, the mathematical assumption underlying HK surrogate modeling is that the bias of the LF simulator has a lower frequency content compared to the response of the HF simulator. Hence, this study proposes to derive the LF simulator as a reduced-order representation of the HF simulator usingdynamicsubstructuring[10]. Inparticular, itisshownhowtoobtainanon-linearreduced-order model as an assembly of reduced-order linear components [28] connected by non-linear springs. 2.Global-sensitivity-based dimensionality reduction. Create a large Experimental Design (ED) fYLF;p;Xg, where YLF;p =fy(1) LF;p;:::;y(NLF) LF;pg,X=fx(1);:::;x(NLF)g,x(n)is a generic sample of the input parameters, and y(n) LF;pis the corresponding LF estimate of the p-quantile of the QoI. Train a non intrusive PCE surrogate model [21], ^MPCE LF;p(x), based onfYLF;p;Xg. Perform a PCE-based GSA [22] w.r.t. the QoI. Select and retain the input random variables that contribute most to the variability (variance) of the QoI and ﬁx the remaining ones to their expected value. In practice, deﬁne X= [X ;x ], whereX are the "important random variables," and x =E[X ]the rest of input parameters. A natural output of GSA are univariate eﬀects of input parameters, which are the expectation of a QoI conditioned to the value of a single input variable (therefore, they provide information about the shape of the surro- gate models). In this context they are used at step 4 for deﬁning monotonic fragility models. A detailed summary on PCE-based GSA is reported in Appendix A. 3.Hierarchical Kriging surrogate modeling. Train a second smaller ED fYHF;p;X g, where YHF;p = fy(1) HF;p;:::;y(NHF) HF;pg,X =fx(1) ;:::;x(NHF) g,x(n) is a generic sample of the selected model parameters, andy(n) HF;pisthecorrespondingHFestimateofthe p-quantileoftheQoI.TrainaHKsurrogate ^MK HF;p(x ) based onfYHF;p;X gusing the PCE surrogate model of the LF structural simulator as trend function [29, 23]. A short summary of HK surrogate modeling is reported in Appendix B. 4 DRAFT4.Probabilistic characterization of the QoI. Compute the statistics of the QoI of interest (including fragility analysis) via Monte-Carlo-based UQ analysis using the HK surrogate model. Speciﬁcally, given an ordered sequence of quantiles P= [p1;:::;pn], steps 2-3 can be implemented sequentially to obtain n maps (i.e., surrogate models) between YPandX. Therefore, this nordered sequence of surrogate models provides the complete discretized CDF of the QoI w.r.t. the input X=x, i.e.F(yPjx). Observe that such description is suﬃcient for a fast inversion of Eq.6. Finally, fragility models of the QoI w.r.t. the importantrandomvariables( X )arecomputedbyMCS.Observe, thatnotallthe marginal fragilityw.r.t. the single important variables are necessarily monotonous. Therefore, following the current practice, only the monotonic increasing fragility are suitable for the PEER-PBEE framework. The univariate eﬀects of QoI quantiles computed at step 2 guides this ﬁnal selection. Despite the use of stochastic simulators (Eq.3), the outlined framework has the critical advantage of providing by construction (Eq.6) a—deterministic—map between the input parameters and a given quantile of the dis- tribution of the QoI. The nature of the uncertainty associated with this distribution is aleatory. Speciﬁcally, it is the aleatory component, Z, of the stochastic model, Ms(Zjx). In fact, in the special case of a deterministic simulator any quantile will collapse into the deterministic solution of Mc(x). A cornerstone of the proposed framework is model calibration and validation with experimental data. In fact, the LF simulator MLFcan be characterized by an oversimpliﬁcation of the real problem, phenomenological (and not validated) assumptions, and non-physical parameters. As a result, a rigorous global model validation becomes prohibitive, leading to not negligible model biases. On the other hand, the advantage of LF structural simulators lies in their reduced computational cost, enabling global studies. Few considerations need to be made in the context of earthquake engineering. In this case, an EDPof interest is the QoI; moreover, a vector of IMs is composed by either some of the "important" input parameters (i.e., IM2X ) or by a set of instrumental IMs ( IM = [PGA;PGV;::: ]). In the latter case, a statistical conversion [30] between the ground motion parameters and the instrumental IMs can be easily derived provided a large set of simulations. The proposed framework allows the computation of any statistic of interest for the selected EDP, including classical fragility functions and surfaces. In the fragility analysis context, a key passage of this framework is the GSA (step 2), which provides a rigorous selection of the important input parameters (i.e., IM) w.r.t. to the EDPvariability. Observe, that this condition is explic- itly dependent on the quantile of (major) interest. Therefore, it potentially eliminates the loss of correlation betweenIMandEDPin the tail of fragility functions [31] (if this probability region is of interest). 3 The Rio Torto Bridge case study The Rio Torto Bridge is a RC structure built in the ’70s and located on the motorway between Florence and Bologna, Italy. The bridge comprises two independent roadways, each of 400m length and characterized by a 13-span deck sustained by twelve piers. Three out of thirteen spans accommodate Gerber saddles, which protect the deck from restrained deformation due to pier settlement. Pier heights vary between a minimum of 13:80m (Pier # 12) to a maximum of 41:00m (Pier # 7). Figure 1a provides an overview of Rio Torto Bridge. The Rio Torto Bridge was selected as a case study to investigate the seismic performance of old RC bridges and possible seismic retroﬁtting schemes within the framework of the RETRO project [17]. The seismic response of the bridge was evaluated via HS. In detail, mock-up models of Pier #9 and #11 of 1:2.5 scale were built and tested via HS at ELSA Laboratory of the Joint Research Centre of Ispra, Italy. Figure 1b reports a schematic of the hybrid model of the Rio Torto Bridge. Non-linear springs represent the in-plane stiﬀness of numerical and physical piers activated by the lateral displacement of the deck. A schematic view of tested piers is also reported in Figure 1b. Figure 2a reports a schematic of the experimental setup conceived for Pier #11. In this speciﬁc case, actuators labeled 1A and 3B control the lateral in-plane displacement of the pier whereas, hydraulic jacks labeled 1B and 2B impose vertical forces owing to the deck weight. Figure 2b provides a picture of the entire installation at the JRC of Ispra. The East-West and the North-South components of the Emilia earthquake of May 29, 2012, recorded from Mirandola station and downloaded from the ITACA database [32] were considered as serviceability and ultimate limit state accelerograms (SLS and ULS, respectively). The response spectra of both SLS and ULS accelerograms are shown in Section 5 where the stochastic modeling of the seismic input for the fragility analysis is described. In this study, the uncertainty on the structural simulator parameters is considered negligible w.r.t. the ground motion parameters, therefore, we set Xc=xc∗, andXXa. For a comprehensive description of the HS campaign, the reader is addressed to [18]. ∗This corresponds to assign a delta Probability Density Function, i.e., (xcxc) 5 DRAFT (a) (b) Figure 1: The Rio Torto Bridge case study: a) side view retrieved from [18]; b) hybrid model. (a) (b) Figure 2: HS setup at JRC ispra: a) frontal and top view of Pier #11 retrieved from [19]; b) picture of the setup from top retrieved from [18]. 6 DRAFT (a) (b) Figure 3: HF simulator of the Rio Torto Bridge: a) deck-to-pier connection b) pier discretization and ﬁbre-based cross-sections. 4 Computational modeling of the bridge Two structural simulators of a single lane of the Rio Torto Bridge were implemented, namely, a HF ﬁnite- element and a LF state-space structural simulator. The former was initially implemented in OpenSees [33] to support the design of the experimental campaign. The latter was obtained, instead, via dynamic substructuring of the HF simulator to assemble reduced-order numerical substructures for HS [18]. Accordingly, LF and HF structural simulators show almost the same response in the linear regime. 4.1 High-ﬁdelity ﬁnite-element model The HF simulator was implemented in OpenSees [33]. Linear Bernoulli beam elements were used for the deck, whose response was assumed not to exceed the linear regime. Cylindrical hinges model Gerber saddles allowing forlateralandverticalrotationsbetweendeckelementswhileblockingrelativetorsionalrotations. Displacement Degrees-of-Freedom (DoFs) of both abutments were constrained, but not rotations. A rigid link accounted for the oﬀset distance between the cap beam axis, and the center of gravity of the deck cross-section. Each pier was clamped at the base. Figure 3 provides schematic views of both deck cross-section and ﬁnite-element discretization of piers. All piers were modeled using non-linear ﬁber-based elements that allow for discretizing the longitudinal steel reinforcement of the member cross-section. The zero-tensile-strength Concrete01 OpenSees material was adopted for concrete assuming a compressive strength fpc=11:47MPa and pc=17:16MPa for piers characterized by solid and hollow cross-section columns, respectively. These values were calibrated based on the quasi-static response of Pier #9 and #11 obtained via HS. In both cases, the ultimate concrete compressive strength was computed as fpc;u= 0:88fpcMPa and the Poisson ratio was set to = 0:2. Corresponding strain values were set to pc=0:002andpc;u=0:006, respectively. The Steel02OpenSees material was adopted for steel rebars assuming a tensile strength fy= 360MPa, along with a Young modulus E= 205103MPa and hardening modulus Ep= 0:025E. ThePinching4 OpenSees material was used to model the hysteretic shear response of transverse beam elements. Parameter values were calibrated for a series of cyclic tests conducted on a 1=4scale mock-up model of Pier #12 documented in [34]. The OpenSees section-aggregator object was used to couple bending and shear models within the same force-based beam element. Finally, a proportional Rayleigh damping model was adopted with mass multiplier d= 0:27300and stiﬀness multiplier d= 0:00787. Figure 4 compares the hysteretic loop of the lateral restoring force of Pier #11 obtained via HS (SLS and ULS accelerograms, as builtconﬁguration) to corresponding HF simulator predictions. As can be appreciated, the HF simulator eﬀectively reproduces the lateral restoring force response of Pier #11. A single evaluation of the time history response of the HF simulator subjected to a 15sec accelerogram sampled at 1msec takes about 10min if performed with a standard laptop equipped with an Intel 1:80GHz i7-8565U CPU and 16GB RAM or similar. 4.2 Low-ﬁdelity state-space model A linearized version of the HF simulator supported the computation of reduced-order mass and stiﬀness ma- trices of deck and piers of the LF simulator via component-mode synthesis [28], which is a speciﬁc dynamic 7 DRAFT (a) (b) Figure 4: Measured and simulated hysteresis loops of the lateral restoring force of Pier #11: a) SLS; b) ULS. substructuring method [10]. In detail, the reduced-order model of the deck was obtained via static condensation retaining only lateral displacement DoFs of piers and Gerber saddles, ( _u=v _v=m1(lrcvkutag(t));(8) whereuandvare the 18 retained displacement and velocity DoFs, respectively; k,candmare reduced-order stiﬀness, damping and mass matrices, respectively, whereas tis the seismic mass vector. The latter multiplies the seismic accelerogram ag(t)to obtain the seismic loading. The 1812Boolean matrix lcollocates the 12 pier restoring forces stored in r=fr1;:::;r 12gto the related state-space equation. A non-linear single-DoF system based on a Bouc-Wen restoring force was calibrated to mimic the response of the corresponding pier in the HS structural simulator,8 >< >:_ui=vi _vi=m1 i[gi(t)civiritiag(t)] _ri= [Ai+ ( isign(rivi) i)jrijn]vi;(9) whereui,vi, andriare the three state variables, namely, displacement, velocity, and hysteretic restoring force, whereasgi(t)represent the deck reaction force of the i-th pier. Mass miand seismic mass tiwere obtained via staticcondensationaswell. Foreachpier, thetoplateraldisplacementDoFwasretainedasmasterDoF,whereas the others were condensed. Consistently with the HF simulator, viscous damping ciwas computed assuming a proportional Rayleigh damping with mass multiplier d= 0:27300and stiﬀness multiplier d= 0:00787. The parameters A, , andnrefer to the Bouc-Wen model [35], which describes the evolution of the hysteretic restoring force ri. As analogously done for damping, mass, and seismic mass, the parameter A, which represents the initial stiﬀness of the non-linear spring, was obtained via static condensation of the HF pier model. On the other hand, nwas assumed constant and equal to one, whereas and were calibrated to match the static response of the corresponding OpenSees pier model subjected to ten sinusoidal displacement cycles producing 0:5% drift. Regarding Pier #11, Figure 5 compares restoring force, dissipated energy history and hysteresis loop of the calibrated Bouc-Wen spring to the corresponding OpenSees simulation. As can be appreciated from the ﬁgure, the Bouc-Wen spring fairly reproduces both restoring force peak and dissipated energy of the corresponding OpenSees model. Figure 6 compares the dynamic response of Pier #11 obtained from time history analyses of LF and HF structural simulators considering the ULS accelerogram of the HS campaign described in Section 3. A single evaluation of the time history response of the LF simulator subjected to a 15sec accelerogram sampled at 1msec performed with a standard laptop equipped with an Intel 1:80GHz i7-8565U CPU and 16GB RAM takes about 20sec. 5 Stochastic modeling of the seismic input The lateral motion of the Rio Torto Bridge is characterized by a relatively long period of vibration. Conse- quently, its reliability is not strongly aﬀected by classic broadband motions. On the contrary, it is vulnerable to pulse like motions with frequencies close to the structural ones. Generally, the presence of a strong pulse component is more likely in near ﬁeld motions. As a consequence, this study employs the stochastic model proposed by Dabaghi and Der Kiureghian [20], which is designed for such type of excitations. 8 DRAFT (a) (b) (c) Figure 5: LF and HF predictions of the cyclic response of Pier #11: a) restoring force, b) dissipated energy, c) hysteresis loop. (a) (b) (c) Figure 6: LF and HF simulations of the seismic response of Pier #11: a) displacement, b) dissipated energy, c) hysteresis loop. The model is composed by a residual broad-band process (as done in [7]), and a modiﬁed version of the Mavroeidis, pulse model with random parameters [36], Vpul(t) =Mpul(t;Xpul) =1 2Vpcos 2tTmax;p Tp + Dr Tp 1+cos2 tTmax;p Tp ,ti;p<ttf;p: (10) whereXpul= [Vp;Tp;Tmax;p; ;]. Speciﬁcally, VpandTpare pulse amplitude and period, respectively, represents the number of oscillations within the pulse, is the phase angle shift w.r.t.the time modulation function, and Tmax;pis a location parameter for the peak of the excitation. Drhas been introduced in [20] to guarantee zero residual displacement at the end of the excitation, and it is deﬁned as Dr=VpTpsin (+ )sin ( ) 4(1 2): (11) Finally,ti;p=Tmax;p0:5 Tp, andtf;p=Tmax;p + 0:5 Tprepresent the start and the end of the pulse motion. While the pulse is modeled in the velocity domain, the residual is modelled in the acceleration domain, i.e., ares(t) =a(t)_vpul(t)∗. Since the residual of the acceleration is a broadband motion, the model follows the same principles of [7], i.e., Ares(t) =Mres(t;ZjXres) =q(tjXq)h(tjXh)W(t;Z) h(t;ZjXh); (12) wheredenotes time convolution; Xres= [Xq;Xh];q(tjXq)is a parametric time modulating function with random parameters Xq;h(tjXh)is the impulse-response function of a linear ﬁlter with time varying random parameters, where Xhrepresents a set of time invariant parameters; W(t;Z) =PNt n(ttn)Znis a band- limited white noise process, where Z= [Z1;:::ZNt]is a standard normal Gaussian vector; and h(t;ZjXh)is the variance of the convoluted process (i.e., the numerator of Eq.12). The modulating function proposed in [20] for near-ﬁeld motions is q(tjXq) =( C t Tmax;q ,0<ttmax;q Cexp [ (tTmax;q )],tmax;q<t;(13) ∗lower case letter are used to indicate the single excitation; moreover, notice the time derivative on the pulse component. 9 DRAFTTable 1: Selected events associated with the Emilia earthquake of May 29, 2012 (Italy). Event ID Time [h:m:s] MWD[km]Station code Latitude [deg] Longitude [deg] IT-2012-0011 07:00:02 64.1 MRN 44.878231 11.061743 IT-2012-0011 07:00:02 64.1 MIRE 44.878212 11.061747 IT-2012-0011 07:00:02 69.3 T0814 44.793300 10.969200 IT-2012-0011 07:00:02 64.5 MIRH 44.882400 11.063100 IT-2012-0011 07:00:02 65.1 MIR02 44.886948 11.073198 IT-2012-0011 07:00:02 60.5 MIR01 44.844042 11.071316 IT-2012-0011 07:00:02 611.2 MIR03 44.938400 11.104500 IT-2012-0011 07:00:02 613 MIR04 44.927433 11.178312 IT-2012-0010 10:55:56 5.56.8 T0819 44.887300 10.898700 IT-2012-0011 07:00:02 611.3 T0813 44.877800 11.199200 IT-2012-0011 07:00:02 614.4 T0800 44.848600 11.247900 IT-2012-0011 07:00:02 66.1 SAN0 44.838000 11.143000 IT-2012-0011 07:00:02 68.6 MIR08 44.916900 11.089500 IT-2012-0011 07:00:02 69.9 T0802 44.875000 11.181600 IT-2012-0011 07:00:02 610.7 T0818 44.934800 11.030400 IT-2012-0011 07:00:02 614.3 T0811 44.783800 11.226500 whereXq= [ ; ;C;T maxq]are the parameters of the function. Speciﬁcally, and are shape parameters controlling the ramping and decreasing phase of the residual, Cis a scale factor, and Tmax;qis a location parameter deﬁning the peak of the residual excitation. Moreover, it can be shown that a one-to-one mapping exists between these parameters and D595;Ia;res=T95;Ia;resT5;Ia;res,T30;Ia;res, andIa;res—whereTp;Ia;resis the time corresponding to the p% of the cumulative Arias intensity of the residual, Ia;res. It follows that Xq can be written as Xq= [Ia;res;T30;Ia;res;D595;Ia;res]. The ﬁlter is deﬁned by h(tjXh) =!f(t)q 12 fexp [!f(t)ft] sinh !f(t)q 12 fti : (14) In Eq.14, the main frequency, !f(t), is evolving with time linearly. Speciﬁcally, !f(t) =!mid+ _!(tT30;Ia), where!midis the frequency at time T30;Iares(i.e., at 30 % of the cumulative Arias intensity of the residual) and _!is the rate of change of the frequency with time. The damping of the ﬁlter, f, is considered time invariant. It follows thatXh= [!mid;_!;]. For further details the reader should refer to [20]. Finally, the artiﬁcial ground motion model can be written as A(t) =Ma(t;ZjXa) =Mres(t;ZjXres) + _Mpul(t;Xpul) =Ares(t) +_Vpul(t); (15) whereXa= [Xres;Xpul]. In order to ﬁnd parameter ranges consistent with the seismic hazard characteristic of the site, the artiﬁcial ground motion model was calibrated against real seismic records following the procedure explained in [20]. In detail, the seismic records of real events were selected based on the disaggregation analysis of the PGA with probability of exceedance of 2% in 50years [37]. Then, it was found that only pulse-like ground motions were likely to cause damage to the structure. Consequently, only pulse-like motions were retained and used for the calibration of the artiﬁcial ground motion model. It follows that the analysis and the result of this case study focus on pulse-like excitations. Table 1 reports the list of the 16selected events, whose corresponding accelerograms were downloaded from the ITACA database [32]. All events belong to the Emilia earthquake occurred on May 29, 2012 (Italy). Each of the 16events listed in Table 1 comprise N-S and E-W acceleration records, which were decomposed into principal components as explained in [38]. Velocity pulses where found on12over 162components only, which were retained for the calibration. Figure 7 reports displacement, velocity and acceleration response spectra of the 12retained components as well as the response spectra of the SLS and ULS accelerograms used for the RETRO experimental campaign. The response spectra of the selected records, which correspond to the gray lines, fall within the range established by the response spectra of SLS and ULS accelerograms adopted during experiments and, therefore, utilized to calibrate the HF simulator of the Rio Torto Bridge. The artiﬁcial ground motion model with pulse-like motion was calibrated for each of the 12retained records following the procedure described in [20]. The parameter estimates are reported in Table 2. 10 DRAFT (a) (b) (c) Figure 7: Selected records with pulse-like components: a) acceleration; b) velocity; c) displacement response spectra. Greylinesrefertosingleretainedrecordswithpulse-likecomponentwhereasblackdashedlinesindicate the corresponding 95% conﬁdence interval. Blue and red lines indicate SLS and ULS records used for the HS campaign, respectively. The periods of the ﬁrst ﬁve eigenmodes of the Rio Torto Bridge estimated with the HF simulator are also reported. Table 2: Parameters of Dabaghi’s model identiﬁed on the pulse-like motions. minandmaxbounds, average values, standard deviations and coeﬃcients of variations CoVare computed over the 12selected records. RecordVpTp tmax;pIa;resT30;Ia;resD595;Ia;res!mid _! 10.2822.0941.2722.9062.7391.029 2.040 7.095 29.327 0.8230.338 20.2822.1973.2622.3482.1330.777 2.260 7.490 33.884 0.9500.444 30.2922.0721.3022.8932.5460.932 1.975 7.020 25.006 0.4930.304 40.2762.1393.3272.4341.9470.619 2.105 7.115 23.464 1.2640.405 50.3032.3341.4092.8312.6150.515 2.070 7.040 20.953 0.4720.423 60.2492.3942.9852.4911.8340.297 2.055 7.465 19.566 1.3680.579 70.3632.6023.5072.4492.5800.897 2.270 7.265 39.586 0.1350.585 80.1962.4251.6872.9173.2510.716 2.385 7.890 43.533 -0.3440.462 90.3922.2824.5302.3502.8620.881 2.305 7.185 32.215 0.2120.306 100.2921.9492.6162.3391.4720.408 1.530 4.195 23.015 1.9480.408 110.2181.7623.0692.4181.7480.335 1.390 3.115 33.868 0.9490.378 120.2052.6881.9242.2123.3520.298 2.025 6.720 30.052 -0.8750.349 min0.1961.7621.2722.2121.4720.297 1.390 3.115 19.566 -0.8750.304 max0.3922.6884.5302.9173.3521.029 2.385 7.890 43.533 1.9480.585 0.2792.2452.5742.5492.4230.642 2.034 6.633 29.539 0.6160.415 0.0560.2541.0010.2490.5730.255 0.286 1.378 7.159 0.7420.089 CoV0.2010.1130.3890.0980.2360.398 0.141 0.208 0.242 1.2030.214 Figure 8 compares 50realizations of artiﬁcial records to the reference real record (with parameters listed in the second entry of Table 2). 6 Seismic fragility analysis of the bridge Consistently with the notation introduced in the previous sections, input variables and QoIs of for the fragility analysis read, YLF;p=fup;1;:::;up;12g,XXa=fVp;Tp; ;;Tmax;p;Ia;res;T30;Ia;res;D595;Ia;res;!mid;g(16) whereup;irepresents the p-quantile of the peak of the lateral drift of pier i-th. A GSA of pier lateral drift was performed to retain only the relevant parameters of the artiﬁcial ground motion model. In the proposed example, the evaluation cost of a single run of the HF structural simulator is about 10min, that is, about 30 times larger than the evaluation cost of the LF simulator. Moreover, the evaluation of a QoI quantile relies on the evaluation of the structural simulator response for 100realizations of baseline noise of the artiﬁcial ground motion model. Accordingly, PCE-based Sobol’ sensitivity indices were computed for the LF simulator of the bridge. An initial ED of 200samples was evaluated considering pier drift quantiles P=f0:05 : 0:05 : 0:95g, 11 DRAFT (a) (b) (c) Figure 8: Response spectra of 50realizations of the artiﬁcial ground motion model calibrated against the second record of Table 2: a) acceleration; b) velocity; c) displacement response spectra. Red lines refer to the reference record, grey lines represent the 50realizations of the artiﬁcial ground motion model, and black dashed lines indicate the corresponding 95% conﬁdence interval. Table 3: First-order Sobol’ indices of lateral drift peak of Pier #11 (QoI). PlooVpTp tmax;pIa;resT30;Ia;resD595;Ia;res!mid 0.1000.0260.5060.4160.0140.0000.0000.000 0.002 0.000 0.0010.000 0.2000.0110.5140.4120.0190.0000.0000.002 0.001 0.000 0.0040.001 0.3000.0120.5200.4100.0170.0010.0000.005 0.000 0.000 0.0070.002 0.4000.0100.5340.3810.0160.0000.0000.016 0.000 0.000 0.0130.004 0.5000.0120.5370.3620.0160.0000.0010.026 0.000 0.000 0.0170.009 0.6000.0120.5450.3320.0140.0000.0010.040 0.000 0.001 0.0230.012 0.7000.0120.5470.2980.0140.0000.0000.055 0.000 0.003 0.0320.018 0.8000.0120.5390.2650.0130.0000.0000.078 0.000 0.006 0.0390.023 0.9000.0260.5200.2140.0160.0010.0000.115 0.000 0.010 0.0520.040 according to step 2 of the proposed framework. Such task entailed 20;000evaluations of the LF simulator, that is, about 6;667min of computational cost. Since GSA is used to uncover the inner working of the simulator, in this stage, it is legitimate to assume that all input variables are independently distributed. Tables 3 and 4 report PCE-based ﬁrst-order and total Sobol’ sensitivity indices related to the lateral drift of Pier #11 for the quantilesP=f0:1 : 0:1 : 0:9g. Leave-one-out errors loo, which are also reported, highlight fairly accurate PCEs. Similar errors were obtained up to 0:99quantiles for all piers. It is noteworthy that, for regression- based PCE, the calculation of leave-one-out errors requires the evaluation of a single surrogate model based on the entire ED [39]. The contribution of Vp,Tp,Ia;res,!midandto the variability of drift quantiles was systematically dominant for all piers. Constant average values reported in Table 2 were considered for all other parameters in the following calculations. Among retained parameters, univariate eﬀects highlighted monotonic increasing relationships with the QoI for Vp,Ia;resandonly. The results of Figure 9 refer to Pier #11 but the same trend was systematically observed for all piers. Accordingly, VpandIa;reswere employed as IMs for computing the fragility analysis. Following step 3 of the framework, a reﬁned ED of 10 samples was evaluated using the HF simulator of the bridge considering, X =fVp;Tp;Ia;res;!mid;g, (17) x =E f ;;Tmax;p;T30;Ia;res;D595;Ia;res;!midg : (18) 12 DRAFTTable 4: Total Sobol’ indices of lateral drift peak of Pier #11 (QoI). PlooVpTp tmax;pIa;resT30;Ia;resD595;Ia;res!mid 0.1000.0260.5390.4710.0240.0020.0010.012 0.003 0.000 0.0050.007 0.2000.0110.5360.4480.0310.0020.0030.010 0.003 0.003 0.0110.006 0.3000.0120.5360.4350.0280.0030.0030.012 0.002 0.005 0.0150.009 0.4000.0100.5460.4080.0260.0010.0020.024 0.001 0.002 0.0210.008 0.5000.0120.5490.3880.0270.0000.0030.033 0.000 0.000 0.0210.012 0.6000.0120.5550.3580.0230.0020.0020.048 0.001 0.002 0.0290.017 0.7000.0120.5540.3230.0250.0010.0030.062 0.002 0.004 0.0400.022 0.8000.0120.5470.2870.0220.0010.0050.090 0.001 0.011 0.0480.027 0.9000.0260.5260.2300.0250.0050.0040.124 0.003 0.012 0.0600.046 (a) (b) (c) Figure 9: Univariate eﬀects for the lateral drift peak of Pier #11 (QoI) for the quantiles reported in Tables 3 and 4. Lighter curves refer to quantile p= 0:1while darker curves refer to quantile p= 0:9. For each sample of the reﬁned ED, 100baseline noise realizations were considered entailing a total of 1;000 evaluations of the HF structural simulator, that is, about 10;000min of computational time. Then, HK surrogates were trained for diﬀerent ED sizes. Accordingly, HF- nindicates a HK surrogate computed with nsamples of the HF structural simulator response quantiles. Following step 3 of the framework, the PCE surrogate trained on the LF simulator was used as a trend function for the HK surrogates. A Matérn 5=2 correlation function R(XX0j)and second-order polynomial trend function (X)were adopted. According to step 4 of the framework, the fragility models reported in Figures 10 and 11 were computed via Monte-Carlo-basedUQforwardanalysisusingsurrogatemodelsandanEDPdriftthresholdvalue y= 0:01. This value corresponds to a medium damage state characterized by the onset of concrete spalling in the transverse beams, as demonstrated in [34, 40, 41] by means of numerical simulations and experimental tests. Speciﬁcally, the fragility curves reported in Figure 10 were obtained from the PCE surrogate of the LF simulator and three HK surrogates of the HF simulator. In the latter case, the three HK surrogates were obtained considering 8, 9, and 10 samples of the HF simulator response; hence, they are indicated as HF-8, HF-9, and HF-10, respectively. The red shaded area indicates the 95% conﬁdence interval of the fragility curve computed with the HF-10 HK surrogate. As can be observed, fragility curves approach a stable shape already with 8samples of the HF simulator response. The motivation behind the selection of Vpas IM for the fragility curves is twofold. First, Vpis associated with the larger Sobol’ sensitivity index and, therefore, it is the most sensitive parameter of the LF simulator (see Tables 3 and 4 in this regard). Moreover, according to Figure 9a, Vpshows a monotonic increasing relationship with the selected EDP. On the other hand, the fragility surface of Figure 11 was obtained for the HK surrogate indicated as HF-10. Owing to their substantial contribution to the variability of the EDP demonstrated by Sobol’ indices, VpandIa;reswere selected as IMs. Moreover, according to Figure 9, both these two parameters show a monotonic increasing relationship with the EDP. In order to demonstrate the accuracy of the proposed framework, both fragility models were compared to crude Monte Carlo estimates of the failure probability of the HF simulator computed for two validation points. In this regard, the red diamonds of Figures 10 and 11 indicate mean value and 95% conﬁdence interval of failure probability estimates obtained with 100samples∗of the HF simulator response—i.e., circa 1;000min per point of computational time. In detail, the validation point for the fragility curve (Figure 10) assumes Vp= 0:3while the validation point for the ∗Notice that with 100 simulations the mean estimate is stable, while the conﬁdence bounds should be considered as indicative 13 DRAFT Figure 10: Fragility curves for the lateral drift peak of Pier #11 (QoI). (a) (b) (c) Figure 11: Fragility surface for the lateral drift peak of Pier #11 (QoI): a) fragility surface; b) section along Vp; c) section along Ia;res. fragility surface (Figure 11) assumes Vp= 0:3andIa;res= 0:6. For the sake of clarity, two orthogonal sections of the fragility surface intersecting the validation point are also included in Figure 11. For both fragility models, Monte Carlo estimates are consistent with fragility model predictions. Although HK surrogates were obtained assuming independent uniformly distributed ground motion parame- ters, fragility models can be easily re-sampled considering a joint PDF calibrated for the speciﬁc seismic hazard. Also, provided with such simulation framework, ground motion parameters can be mapped to instrumental IM as illustrated in [30]. 7 Conclusions This study introduced a computational framework for eﬃcient and accurate seismic fragility analysis based on a combination of artiﬁcial ground motions and surrogate modeling of multiple stochastic simulators characterized by diﬀerent levels of ﬁdelity. The eﬀectiveness of the proposed framework is demonstrated throughout a bench- mark case study, which consists of a reinforced concrete bridge vulnerable to ground motions with pulse-like components. Speciﬁcally, surrogate modeling is applied at the level of quantiles of a scalar response quantity of interest, which is evaluated either with an expensive-to-evaluate high-ﬁdelity or a cheaper low-ﬁdelity stochastic simulator. As a result, two fragility models associated with a medium damage state are computed for the lateral drift of one pier of the bridge based on 1;000evaluations of the high-ﬁdelity simulator plus 20;000evaluations of the low-ﬁdelity simulator. Accuracy of failure probability estimates is demonstrated by comparison with Monte Carlo simulations of the high-ﬁdelity simulator. The main ﬁndings of the benchmark study, which can be reasonably generalized to seismic fragility analysis of bridges, are that: •100artiﬁcial ground motions are suﬃcient to attain fairly accurate quantile estimates in the range 0:05 0:95for a given sample of the stochastic simulator parameters (both low- and high-ﬁdelity). •Polynomial-chaos-based global sensitivity analysis facilitates the selection of the important parameters of the artiﬁcial ground motion model to be used as intensity measures for the fragility analysis. •Dynamic substructuring and hierarchical kriging provide a mathematically sound framework for balancing computational cost and accuracy of structural analysis. 14 DRAFTIt is critical to remark that the domain of validation of the high-ﬁdelity simulator deﬁnes the perimeter of any validated output of the framework. An important future research direction that directly builds upon this framework is the optimal design of experimental campaigns to meet a target degree of accuracy of fragility models. Appendix A. PCE-based Sobol’ sensitivity indices Assume that the vector of input variables Xhas supportDXand follows an independent joint PDF fX(X) =QM mfXm(xm)wherefXmis the marginal PDF of the m-th input variable. Any square-integrable mapping Y=M(X)w.r.t. the probability measure associated with fX, can be written as a sum of functions of increasing dimension as [42]: M(X) =M0+MX m=1Mm(Xm) +X 1mlMMm;l(Xm;Xl) +:::+M1;2;:::;M (X); (A.1) or equivalently: M(X) =M0+X u6=;Mu(Xu); (A.2) whereM0is the mean value of Y,u=fm1;:::;msgf 1;:::;Mgare index sets, and Xudenotes a subvector ofXcontaining only the components indexed by u. The number of summands in the above equation is 2M1. The uniqueness and orthogonality of Sobol’-Hoeﬀding decomposition allow for the following decomposition of the variance DofY: D=Var[M(X)] =X u6=;Du (A.3) whereDudenotes the partial variance: Du=Var[Mu(Xu)] =E M2 u(Xu) : (A.4) The Sobol’ index Sucan be deﬁned as the fraction of the total variance Duthat corresponds to the set of input variables indexed by u: Su=Du D: (A.5) Byconstruction,P u6=;Su= 1. First-orderindices S(1) mdescribetheinﬂuenceofeachparameters Xmconsidered separately. Second-order indices S(2) mldescribe the inﬂuence from pairs of parameters fXm;Xlgnot already accounted for by XmorXlseparately. High-order indices describe combine inﬂuences from larger sets of parameters. The total sensitivity indices S(tot) mrepresent the total eﬀect of an input variable Xmaccounting for its main eﬀect and all interactions with other input variables. It follows that S(tot) m = 1Sm, where Smis the sum of all Suwithunot including m. Sobol’ indices can be evaluated by Monte Carlo simulation [43], which requires O(103)model evaluations for each index Su. Sobol’ indices can be obtained analytically at no additional cost than computing the PCE. A concise description of PCE-based Sobol’ sensitivity indices estimation is given herein; for further details, the reader is referred to [44, 45]. PCE relies on the decomposition of Y=M(X), as a linear superposition of non-linear functions as follows: bY=MPCE(X) =X 2Ay (X); (A.6) wheref ; 2Agis a set of multivariate polynomials that are orthogonal w.r.t. the input vector with independent components XfX(X) =QM m=1fXm(Xm), = ( 1;::: M)is a multi-index that identiﬁes the polynomial degree in each of the input variables, and y denotes the corresponding polynomial coeﬃcient (coordinate). Following the orthonormality condition, E MPCE(X) =y0. For practical purposes, the inﬁ- nite sum in Eq.A.6 needs to be truncated to a ﬁnite series. This is commonly achieved by maximum-degree or hyperbolic norm truncation (for more details, see [21]). In order to compute the PCE, ﬁrst, a so-called ED, consisting of a set of realization of the input vector X=fx(1);:::;x(N)gand the corresponding model 15 DRAFTevaluationsY=fy(1);:::;y(N)g, is generated. Then, the set of coeﬃcients y is estimated by minimizing the expected mean-square approximation error on the ED by solving: by = arg min y E YbY2 (A.7) In the present application, the ED was formed by sampling the input variable space with a Sobol’ low- discrepancy sequence. The minimization in Eq.A.7 was solved using the hybrid least angle regression method originally proposed in [21]. It is straightforward to obtain the Sobol’ decomposition of Yin an analytical form by observing that the summands MPCE u(Xu)in Eq.A.2 can be written as: MPCE u(Xu) =X 2Auy (Xu); (A.8) whereAu=f 2A : k6= 0if and only if k2ugdenotes the set of multi-indices such that [Au=A. Consequently, due to the uniqueness of Sobol’-Hoeﬀding decomposition, there is an analytical expression of MPCE uof Eq.A.8, which serves as a proxy of Mu. The analytical expression of the total variance of a PCE is given by [46]: D=Var MPCE(X) =X 2Ay2 : (A.9) Similarly, the partial variance Dureads: Du=Var MPCE u(Xu) =X 2Auy2 : (A.10) Accordingly, the Sobol’ indices of any order can be approximated by a simple combination of the squares of the PCE coeﬃcients by substituting Eqs.A.10 and A.9 in Eq.A.5. For instance, the ﬁrst-order Sobol’ indices, which describe the inﬂuence of each input variable Xmconsidered separately, read: S(1) m=P 2Amy2 P 2Ay2 ,Am=f 2A: m>0, m6=l= 0g (A.11) whereas the total Sobol’ indices, which represent the total eﬀect of an input variable Xmaccounting for its main eﬀect and all interaction with other input variables, are given by: S(tot) m=P 2Atotmy2 P 2Ay2 ,Atot m=f 2A: m>0g: (A.12) Sobol’ indices provide quantitative insight on the importance of an input variable. However, they do not include information about the direction in which an input variable aﬀects the model response Y. So-called univariate eﬀects can answer this question [47, 48]. A univariate eﬀect is the expectation of Yconditioned to the value of a single input variable. Univariate eﬀects have a closed analytical form for PCE models, closely related to the ﬁrst-order Sobol’ decomposition, M(1) m(Xm) =X 2Amy (Xm),Am=f 2A: m>0; m6=l= 0g (A.13) In this study, computations of both PCE and Sobol’ indices were performed using UQLab, which is a MATLAB toolboxforUQdevelopedbytheChairofRisk, SafetyandUncertaintyQuantiﬁcationofETHZurich[49,50,51]. Appendix B. HK surrogate modeling LetfMc; c= 1;:::Cgbeing a series of simulators sorted by increasing level of ﬁdelity, and ycthe associated output,yc=Mc(x). In addition, let ^MK cdeﬁne the Kriging surrogate model of the simulator Mc(x). A Kriging surrogate model is deﬁned as an inﬁnite collection of jointly Gaussian random variables (i.e., a Gaussian stochastic process). Any ﬁnite sample of the process is a Gaussian random vector, which is completely 16 DRAFTdeﬁned by the multivariate Gaussian distribution [52]. In this context, the input Xis sampled at Ncdistinct locations within the support DX, and the corresponding scalar output aren y(n) c=Mc x(n) c ;n= 1;:::;Nco . TheNcdistinct input samples are collected in the set Xc= X(1);:::;X(Nc) Tand the output in the set Yc=n y(1) c;:::;y(Nc) coT . Here, the set Ycis considered as realization of a Gaussian random vector, Yc=  Y(1);:::;Y(Nc)T, with mean c(X), and covariance matrix c(X). Therefore, the kriging surrogate model of thecsimulator can be deﬁned as [53, 54, 52] Mc(X)^MK c(X) =c(X) +2 cZ(X); =fc(X)T c+2 cZ(X);(B.1) where cis a vector of regression coeﬃcients, fc(X)is a vector collecting a series of basis functions, 2 cis the variance of the process, and Z(X)is a zero-mean, unit-variance stationary Gaussian process. Therefore, Z(X)is fully determined by the correlation function R(X;X0j) =R(XX0j)between two distinct points (X;X0)in the input space, where is a set of hyper-parameters. Given the input ED Xcand the output set Yc, c;c, andcare generally determined by generalized least squares as follow [54, 55]: ^c= arg minAD1 2log (det (R)) +n 2log 22 c +n 2 ; (B.2) ^ (^c) = FT cR1Fc1FT cR1Yc; (B.3) ^2 c(^c) =(YcFc^ c)TR1(YcFc^ c) Nc; (B.4) whereRn;m=R(Xn;Xm)isthecorrelationmatrixbetweenthe Ncsamples,and Fc= [fc(X(1));:::;fc(X(Nc))]T is the information matrix. Next, given a desired test point x, the random variable of the unobserved output Y cjXc;Yc=MK c(xjXc;Yc)has conditional Gaussian distribution [54] with parameters yc=Eh Y cjXc;Yci ; =fc(x)T^ c+rT(x)R1 YcFc^ c ; (B.5) 2 yc= Varh Y cjXc;Yci ; = ^2 ch 1rT(x)R1r(x) + FT cR1r(x)fc(x)T FT cR1Fc1 FT cR1r(x)fc(x)i ; (B.6) wherer(x) =^c) = [R(xx(1)j^c);:::;R (xx(Nc)j^c)]Tis the cross-correlation vector between the point xand each of the points of Xc. In this study, we use a nested combination of surrogate models to deﬁne a HK surrogate [23]. Speciﬁcally, a PCE surrogate of the LF simulator is used as a trend for the HF simulator. For the two levels c2[LF;HF], the HK predictor at an unobserved point xcan be written as [29]: y HF=^MPCE LF(x)^ HF+r(x)TR1(YHFFLF^ HF); (B.7) where ^MPCE LF(x)is PCE surrogate of the LFsimulator, YHFis the vector of output from the HFsimulator, FLF= [MLF(x1);:::;MLF(xNc)], and ^ HFis a constant. The variance of the HK predictor is similar to Eq.B.6, i.e. 2 y HF= ^2 HF 1r(x)TR1r(x) + FT LFR1r(x)^MPCE LF(x)T FT LFR1FLF1 FT LFR1r(x)^MPCE LF(x) : (B.8) In this study, computations of kriging surrogates were performed using UQLab, which is a MATLAB toolbox for UQ developed by the Chair of Risk, Safety and Uncertainty Quantiﬁcation of ETH Zurich [49, 56]. Acknowledgments The authors wish to acknowledge the Chair of Structural Dynamics and Earthquake Engineering (Prof. Dr. Božidar Stojadinović) and the Chair of Risk, Safety and UQ (Prof. Dr. Bruno Sudret) of ETH Zurich, Switzerland, for having supported this research. 17 DRAFT