Evaluation of the variability contribution due to epistemic uncertainty on constitutive models in the definition of fragility curves of RC frames

Abstract In the framework of uncertainty propagation in seismic analyses, most of the research efforts were devoted to quantifying and reducing uncertainties related to seismic input. However, also uncertainties associated to the definition of constitutive models must be taken into account, in order to have a reliable estimate of the total uncertainty in structural response. The present paper, by means of incremental dynamic analyses on reinforced concrete frames, evaluates the effect of the epistemic uncertainty for plastic-hinges hysteretic models selection. Eleven different hysteretic models, identified based on literature data, were used and seismic fragility curves were obtained for three different levels of maximum interstorey drift ratio. Finally, by means of analysis of variance techniques, the paper shows that the uncertainty associated to the hysteretic model definition has a magnitude similar to that due to record-to-record variability.


Introduction
Performance-Based Earthquake Engineering (PBEE) probabilistic procedures, aim at providing a framework for the systematic treatment of all the uncertainties involved in seismic engineering problems [1,2]. These uncertainties are mainly related to seismic actions, typically named record-to-record variability, and nonlinear structural behaviour in dynamic conditions. The latter uncertainties can be associated to different factors: (a) definition of material properties; (b) material hysteretic behaviour; (c) definition of backbone curves and cyclic behaviour of sections; (d) definition of numerical models (geometry uncertainty, load position, mass distribution, damping, description of damage mechanisms that can lead to collapse, etc.). These uncertainties are sometimes grouped under the label of modelling uncertainty [3].
Most of the vast literature on the topic concerns the quantification of the uncertainty in structural response due to record-to-record variability [4][5][6][7][8][9][10][11][12][13][14]. Whereas the level of knowledge on the effects of modelling uncertainty is much more limited. For instance, for many years, in seismic analyses, properties of materials, constitutive and geometric models, were assumed as deterministic [15][16][17][18]. Recently, some researchers investigated the effects of uncertainties in the definition of materials properties [19][20][21] or constitutive models [3,22], but the consequences of epistemic uncertainties on constitutive model selection (i.e. uncertainties related to the assumption of a specific constitutive model in a numerical simulation) were seldom studied.
This fact could lead to the false conclusion that effects of the latter are less important, or even negligible, than the effects of ground motion variability.
One of the reasons why modelling uncertainty was not investigated as much as record-to-record variability is probably its strong systemdependency (i.e. dependent on constructive technology) which makes not easy to obtain general results. Recently, the importance of considering modelling uncertainties, especially in the definition of collapse capacity has been acknowledged by researchers [23,24]. Vamvatsikos and Fragiadakis [22] evaluated the variability in both seismic demand and capacity, on a nine-storey steel moment resisting frame by means of Incremental Dynamic Analysis (IDA). They adopted moment-rotation relationships with non-deterministic backbones for plastic-hinges in beams, and assessed the seismic performance of their model for several combinations of parameters. The uncertain parameters of the backbones were: the yield moment, the post-yield hardening ratio, the endof-hardening rotation, the slope of the softening branch, the residual moment capacity and the ultimate rotation. The same hysteretic rules were adopted in all models. Vamvatsikos and Fragiadakis have also confirmed that, because of the strongly non-linear nature of the problem, a model with median values of constitutive parameters, does not provide an estimate of the median seismic demand capacity. Dolsek [21] extended the standard IDA procedure by proposing to use a group of structural models, in addition to different ground-motion records, in order to estimate both modelling uncertainty and record-to-record variability. Considering a four-storey RC frame, Dolsek [21] has shown that the influence of modelling uncertainty has a significant influence on seismic collapse capacity. Borgonovo et al. [25] discussed importance measures in seismic probabilistic risk assessment and proposed a method to guide the choices of the modeller to detect critical factors in modelling, with the main goal of producing a suitable-topurpose set of fragility curves obtained starting from a single fragility curve and considering the epistemic uncertainties affecting modelling activities. Celik and Ellingwood [26] by using probabilistic non-linear analyses, have proposed fragilities for Reinforced Concrete (RC) frames designed for gravity loads and concluded that record-to-record variability dominates the total uncertainty in the structural response of lowto-moderate seismic areas. Moreover, they stated that fragilities developed using mean values of structural parameters may be sufficient for earthquake damage and loss estimation in moderate seismic regions. This seems not completely in agreement with the findings in [3]. A similar conclusion is drawn by Kwon and Elnashai [20]. These authors, studying a three story ordinary moment resisting RC frame, indicated that the effect of randomness of material parameters had less importance than the effect of ground-motion variability. Recently, Ugurhan et al. [27], studied the uncertainty in the evaluation of the seismic collapse of modern RC frame structures, concluding that the probability of collapse is sensitive to modelling uncertainties.
The present paper discusses the effects of epistemic uncertainty on the choice of constitutive models for RC frame structures, while the uncertainty in model parameters is not addressed. In particular, in this work, IDA [28] is used in order to estimate the effect of uncertainties in the definition of hysteretic models for plastic hinges in beams and columns of RC frames. By analysing the results of the IDA analyses, we have obtained fragility curves for the maximum interstorey drift ratio, in terms of both Peak Ground Acceleration (PGA) and spectral acceleration for the first natural period of the structures under consideration Sa(T 1 ). Furthermore, by means of ANalysis Of VAriance (ANOVA) techniques, record-to-record variability and constitutive modelling variability are quantified, showing that, at large IDRs, the latter is comparable to the former.

Description of the structure and definition of collapse mechanisms
The structure used as case study is a three-storey three-bay RC frame ( Fig. 1), with the same geometry adopted in [26]. This frame has 30 cm × 30 cm column cross-sections and 23 cm × 46 cm (B × H) beam cross-sections, the length of the bays is 549 cm and the interstorey height is 366 cm for all the three storeys. For the design against vertical loads, a distributed load of 48 kN/m was assumed on the beams. The design for vertical loads complies to Eurocode 2 prescriptions [29].
In order to reproduce the behaviour of different classes of structures, two alternative reinforcement bars (rebars) designs were considered. These designs can be associated to two different types of expected collapse mechanisms (named collapse mechanism CM1 and CM2 in the following). The CM1 design is characterized by weak columnstrong beam failure, that can be representative of existing structures designed without the application of capacity design principles (in particular the strong column-weak beam hierarchy) [30]. The layout of reinforcement bars for this first structure is depicted in Fig. 1b. Beams have 3Ø20 bars at the top and at the bottom of the end sections. In the central portion of the beams there are 2Ø14 and 4Ø16 at the top and at the bottom, respectively. The columns have 3Ø16 rebars for each edge. Based on this design a soft storey failure mechanism is expected. The CM2 reinforcement design is fully compliant to Eurocode 8 prescriptions, thus including capacity design rules. Therefore, in this case a global mechanism of collapse is expected, with plastic hinges at the ends of beams. The beams of this frame contain the same rebars of those of CM1, while total of 12Ø20 are considered for column crosssections (see Fig. 1b). For both reinforcement designs, the shear capacity of columns and beams was assumed higher than the flexural capacity, so to exclude brittle shear failures. Any other possible brittle failure mechanism (e.g. buckling and loss of anchorage of rebars, pullout of the cover, joints failure) was excluded as well.

Finite element modelling
Numerical analyses were carried out using finite element models of the frames. The software OpenSees was employed to this aim [31]. Fig. 1a illustrates the main geometry of the structure considered. Beams and columns were modelled using elastic finite elements with lumped plastic hinges at their ends. At the ground level, columns were fully clamped. Masses corresponding to structural and non-structural dead loads and a fraction of live loads, according to Eurocode 8 criteria, were considered as equivalent distributed masses on the beams. A concrete compressive strength of 36 MPa was assumed as well as an elastic modulus of 30.46 GPa. The adopted steel yielding stress was 500 MPa and its elastic modulus 200 GPa. Both the frames considered had a natural period of T 1 = 1.09 s.
Nonlinearity was introduced in the FE models by considering cyclic nonlinear moment-rotation relationships for the plastic hinges. In order to investigate the effects of the uncertainty in the choice of the hysteretic model, eleven different rules were considered in the study, as discussed in Section 2.3. Furthermore, two types of moment rotation backbone curve were used (see Fig. 2a): (i) bilinear elastic-hardening and (ii) tetra-linear with strength degradation. The bilinear backbone and the first two branches of the tetra-linear model are defined by the same parameters, i.e. yielding moment (M y ) and rotation (θ y ) and peak moment (M peak ) and rotation (θ peak ). The tetra-linear model features a negative stiffness branch after the peak moment capacity. The coordinates of the points of these backbone curves were calculated as recommended by FEMA-273 [32], taking into account the effect of the axial load in static conditions in column cross-sections, and considering values recommended for elements with flexural dominated failure. The residual bending moment (M res ) was assumed as 20% of M peak (M peak in Fig. 2a), as suggested in [32]. The rotation corresponding to the residual moment (θ res ) was considered as a fitting parameter.
The Hysteretic material model available in Opensees was adopted in the analyses. This model describes the moment-rotation relationship based on a tetra-linear backbone curve and five damage parameters related to cyclic degradation (i.e. P X , P Y , D1, D2, β) [33]. The backbone curve of the model (see Fig. 2a) is defined as a function of the coordinates of the end points for the different linear branches, i.e. (M y ,θ y ), (M peak ,θ peak ) and (M res , θ res ). By properly setting the values of these parameters it is possible to represent both the bilinear and tetra-linear backbones discussed above.
As for the damage parameters, P X indicates the pinching factor for rotations (horizontal axis in a moment-rotation diagram) during reloading, while P Y is the pinching factor for moments (vertical axis in a moment-rotation diagram). Unit values for P X and P Y correspond to no pinching while smaller values will introduce the pinching effect. The parameters P X and P Y are typically set together; for RC structural elements with pinching-like degradation mechanisms, the value of the ratio P Y /P X is normally lower than 1.0 and spans from 0.7 (low pinching degradation) to 0.2 (severe pinching degradation). The parameters D1 and D2 are used to model cyclic strength degradation based on ductility and dissipated energy, respectively [33]. In particular, the strength reduction introduced by D1 is proportional to the maximum rotation reached during the loading cycles while the strength reduction associated to D2 is proportional to the hysteretic energy dissipated. The parameter β describes unloading stiffness degradation. Following a ductility-based approach [34], the unloading stiffness K UNL is expressed by the following relationship: where θ max indicates the maximum rotation in the plastic hinge under consideration. The calibration process used to define the values of the aforementioned parameters is discussed in Section 2.4.

Definition of possible hysteretic models
The main scope of the present work is to evaluate and define, from a probabilistic point of view, the effects on seismic fragility curves of the epistemic uncertainty related to hysteretic model selection. Therefore, a set of cyclic constitutive models that represent possible alternatives must be defined. Constitutive models adopted in structural analyses are typically empirically defined based on direct observations or theoretically derived, therefore it is unlikely that a single constitutive model is perfectly representative of the behaviour of structural elements [35]. Even when detailed experimental data is available, the prediction of the seismic structural behaviour remains uncertain as many blind prediction contests have shown [36,37]. Moreover, if the structure under investigation is part of an existing building, uncertainty might be even higher because of limited knowledge on material properties, geometry, detailing, etc. Therefore it is not possible to deterministically define constitutive models for structural members (or cross-sections). This matter is also more complicated when dealing with cyclic (seismic) loading because the expected hysteretic degradation mechanisms must be selected. This might be done using engineering judgment, but, the use of subjective expert opinions might lead to the tendency to include too many models which are often: (i) highly dependent, resulting in redundancy; or (ii) not plausible, resulting in a potentially significant over-estimation of model uncertainty [35,38].
The constitutive models and hysteretic rules adopted here were defined based on literature experimental outcomes. An analysis of the literature on experimental tests of RC elements [39][40][41][42][43][44], allowed to identify their main hysteretic degradation mechanisms. At least five mechanisms, relevant for the definition of the plastic hinge behaviour, were found: two involving stiffness (i.e. unloading stiffness degradation and reloading stiffness degradation), two affecting strength (i.e. cyclic strength degradation and in-cycle strength degradation) and a pinching type degradation (due to slip between concrete and rebars). The pinching-type degradation mechanism has distinctive features and can be considered as independent [45], whereas the other four mechanisms can act in different combinations. Kurtman [46] has identified and quantified the degrading behaviour of 196 RC cantilever elements, using the results of experimental cyclic tests extracted from the PEER Structural Performance Database [39] and from [40]. Furthermore, Kurtman found that the reloading stiffness degradation mechanism was present in almost all the considered cases.
The literature suggests that, among the multilinear models not requiring hysteretic damage parameters (i.e. elastic-plastic model, elastichardening bilinear model, peak oriented model and Clough and Johnston model) the Clough and Johnston model [47] is the most suitable to represent the flexural behaviour of RC structural elements with no or very-limited cyclic degradation. Therefore, Clough and Johnston hysteretic behaviour considering only reloading stiffness degradation was assumed as reference model (labelled as REF code in Table 1) for the calibration of all the other hysteretic models. The Clough and Johnston model was obtained by setting the parameters P X = P Y = 1.0 and D1 = D2 = β = 0.0. Based on the analysis of the literature the eleven different hysteretic models reported in Table 1   Table 2 Values of the damage parameters used for the different constitutive models.  Table 1 Constitutive models adopted in the study with the corresponding damage mechanisms.

Damaging parameters for structures CM1
Constitutive model definition were considered in this study. For each of them Table 1 shows the degradation mechanisms taken into account and the type of backbone curve used. The UNL (unloading) class of models represent the behaviour of elements with reloading and unloading stiffness degrading mechanisms. This class of models were obtained by calibrating the parameter β with the other damaging parameters set as for the REF case. The CYC (cyclic) class has been added with the aim to consider cyclic strength and reloading stiffness degradation mechanisms, identifying a proper value of the parameter D2. The IN-CYC (in-cycle) models reproduce the hysteretic behaviour of degrading hinges for in-cycle strength and reloading stiffness decay mechanisms. This class of models were obtained Fig. 3. Hysteretic moment-rotation relationships used for the various constitutive models for the frame CM1 (a) and CM2 (b). For a better comparison moments and rotations were normalized to their peak values. by using backbone curves with strength degradation (tetra-linear). The rotation θ res was considered as unknown. The pinching-type degradation mechanisms are represented by the class PINCH. In this case the calibration process identified the proper value of the ratio P Y /P X . Finally, a sixth class of models was adopted, in order to consider the contemporary presence of unloading-reloading stiffness degradation and cycle and in-cycle strength degradation (MIX). For the calibration process of this class of models θ res was considered as unknown (see Fig. 2a), having fixed the values of β and D2 on the basis of the experimental evidence collected in [46] and assuming as backbone the curve with strength degradation.

Dissipated energy measure
In order to compare the different hysteretic behaviours under consideration a proper damage measure needs to be defined. Several local damage indices are available in literature [48][49][50][51][52]; some of them depend on the ratio of the total dissipated energy over the total energy available, other parameters depend either on plastic rotation or curvature demand and capacity [53]. Anyway, they seem not fully suitable for the aim of the paper, since they express a ratio between demand and capacity evaluated on a specific constitutive model and therefore do not allow to compare different models. For the aim of the present work a parameter based on the dissipated hysteretic energy was defined. In particular, assuming the Clough and Johnston model as reference, a cumulated Dissipated Energy Ratio (DER) was defined, for a general hysteretic model, as: where E i and E i CJ represent, the hysteretic energies dissipated under a prescribed loading sequence up to the n-th loading cycle, for a general constitutive model and for the Clough and Johnston constitutive model, respectively. Based on experimental data collected in [46], DER values can range from 80 to 90% for structural elements with low to moderate hysteretic damage to a minimum of 55-60% for members with severe hysteretic damage.

Calibration of the hysteretic behaviour of plastic hinges
Moment-rotation backbone curves presented in Sections 2.2 define the behaviour of the plastic hinges under increasing monotonic loading, but they are not sufficient to describe the hysteretic behaviour of RC elements under cyclic reverse loading. To this aim it is necessary to define values for the cyclic damage parameters (see Section 2.3). These were calculated by means of an inverse analysis procedure, based on the results of the wide statistical analysis of the results from [39,40]. This procedure considers a simple cantilever column model, comprising an elastic element and a nonlinear hinge at the base. The hinge has a backbone curve derived from the properties of either the cross-sections of columns, model CM1, or the end sections of beams, model CM2. Cyclic horizontal displacements with increasing amplitude were applied at the top of the cantilever column, for a total of five complete cycles. The cyclic displacement history was set in order to achieve displacement ductility demands of 2, 3, 4, 5 and 6, respectively (see Fig. 2b). The calibration procedure starts from the calculation of the energy dissipated by the cantilever adopting the Clough-Johnston model (REF in Table 1). The parameters of the other ten models were then defined assuming considering target values for the DER: either 0.65 or 0.85. Table 2 reports the values of the damage parameters identified, and then used in the dynamic analyses, for each constitutive model. Table 2 also reports the values of the parameters α 2 and α M (see Fig. 2a where: α 2 = K 2 /K 0 and α M = Mres/Mpeak respectively) derived by the values of moments and rotations. The moment-rotation numerical responses obtained for the various models for the cyclic quasi-static history of displacements are depicted in Fig. 3, overlapped to the response of REF model. For the ease of comparison, moments and rotations are normalized to the values corresponding to the peak-point of the backbone curve. Comparing the hysteretic behaviours depicted in Fig. 3 with the experimental results available in the literature, it is possible to notice that the group of plastic hinges with a DER value of 65% fit very well with the expected behaviour of members with poor detailing and confinement. On the other hand, the hysteretic behaviours calibrated to provide DER = 85% are more representative of elements complying with modern codeprescriptions. The effects of the single degradation mechanisms are clearly evidenced by the hysteretic behaviour of the classes UNL, CYC, IN-CYC and PINCH. The outcomes reported in the following sections will show which degrading mechanisms are able to produce the worst scenario for the investigated models. Based on these findings, it is possible to define a hierarchy of the most severe mechanisms potentially addressing future studies.

Seismic input and dynamic analyses
The dynamic performances of the structures were evaluated by IDAs. These analyses allow to evaluate the response variability caused by constitutive behaviour uncertainty for different levels of displacement demand. The set of thirty ground-motion records selected in [22] were used in the present work. They represent a seismological scenario consisting of a magnitude 6.5-6.9 earthquake at a distance of 15-33 km (fault rupture distance). No record contains near-source effects (i.e. directivity effects). The main characteristics of the records are reported in Table 3, while Fig. 4 shows their elastic pseudo-acceleration response spectra for 5% damping ratio. This set of ground motion records was adopted in order to allow an easy comparison with other literature results. However, it should be noticed that ground-motion records selected according to a seismological scenario [54] (i.e. based on earthquake magnitude and distance only) might lead to an overestimate of record-to-record variability when using IDA analyses, mainly because: i) the reference seismological scenario used for the selection (which can be obtained from hazard disaggregation) might not be representative for all the ground-motion IM values used in IDA; ii) information on spectral shape is not used for the selection [55], therefore the IDA procedure might not be robust against scaling. Other ground-motion selection criteria, such as those proposed by Baker [56] and by Bradley [57], could lead to a more accurate estimate of record to record variability.
In this study two different ground motion intensity measures (IM) were considered: the peak ground acceleration (PGA) and the 5%damped spectral acceleration at the first natural period of vibration, Sa (T 1 ,5%). The PGA is clearly a much less efficient IM that Sa(T 1 ) and therefore will lead to an higher record-to-record variability, nevertheless it was adopted here because it is still often adopted in the literature for the definition of fragility models.
The maximum interstorey drift ratio (IDR) over all the storeys was selected as engineering demand parameter (EDP), in order to define the response of the models. The IDAs were performed by scaling each record starting from a value of IM equal to 0.05 g until structural collapse, which was identified as either the achievement of the collapse IDR value, set equal to 5%, or the onset of dynamic instability. In this latter case, as in [22], flat IDA graph was considered until the collapse IDR (i.e. 5%). P-Δ effects were included in the analyses. A total of 1320 IDA curves were obtained and then analysed. In all the analyses maximum rotations resulted smaller than the ultimate rotation values given by the procedure described by Fajfar et al. [58].

IDA analyses results
The present Section analyses the effects of each parameter and discusses the main outcomes of the uncertainty analysis. In the following, the individual IDA curves for the various ground-motion record are presented for the different types of structural models, together with their 16%, 50%, 84% percentile curves. Fig. 5a reports the results of the IDAs performed on the model CM1 adopting the REF constitutive model by considering PGA as IM, while Fig. 5b shows results in terms of S a (T 1 ,5%). As expected [28,59] the spectral acceleration is more efficient (i.e. provides values less scattered) than PGA for every level of interstorey drift considered. The dispersion in the values obtained using PGA is not negligible also in the quasi-elastic state of the frame. Fig. 5c and d show, for the model CM2 with the REF cyclic behaviour, the PGA and S a (T 1 ,5%) IDA curves, respectively. Comparing the response of structure CM1 (Fig. 5b) and CM2 (Fig. 5d) one can observe that for the REF model, the 16%, 50% and 84% percentile curves have similar values in the IDR range considered.
In order to compare with more detail the variations in the response induced by the different constitutive models, Fig. 6a-d show IDA curves obtained for the frame CM2 for the constitutive models IN-CYC-65, UNL-65, PINCH-65 and MIX-85. In this case there is a relevant difference on the median collapse spectral acceleration (i.e. for IDR = 0.05). It is equal to 1.1 g for the IN-CYC-65, to 0.75 g for UNL-65 and to 0.7 g for the PINCH-65 model. It is worth noting that these three models correspond to the same level of DER = 65%. Furthermore, the model MIX-85, that combines various deterioration mechanisms, provides a lower median collapse spectral acceleration, equal to 0.68 g, even if it has a higher dissipation capacity (DER = 85%) than the other models (DER = 65%). This is not the only case among those considered in the study, therefore we can conclude with the adopted damage index (DER) it might not be possible to define a-priori which is the most  It is important to highlight that the MIX model, produces IDA curves with a markedly irregular trend, most likely because the high complexity of the model leads to a strong sensitivity of the interstorey drift ratio to the ground-acceleration time variations, which cannot be described by a scalar IM. In order to have a comprehensive overview of the main outcomes connected to the adoption of the different constitutive models Fig. 7 shows the median IDA curves for the different constitutive behaviours together with their mean curve, for the different CM. In general curves IDR vs. PGA result more scattered than IDR vs. Sa(T 1 ,5%). The structures have collapse spectral acceleration S a (T 1 ,5%) ranging between 0.35 g and 1.0 g. It is worth noting that the models IN-CYC-65, IN-CYC-85 and MIX-65 in Fig. 7a and b, having a softening behaviour with incyclic strength degradation, produces for structure CM1 very flat curves after IDR around 2% and until collapse (i.e. IDR = 5%). This trends occur because of the arising of a soft-storey mechanism at lowest floor and then the acceleration level activating the maximum capacity of the base plastic hinges is practically the collapse acceleration. Anyway, by comparing the various median IDA curves, is not possible to draw general considerations on the constitutive models producing the dynamic response reaching the highest/lowest structural capacities.

Fragility functions
This section presents fragility curves associated to each constitutive model, for the two different classes of structures (CM1 and CM2) considered in the paper. Following a consolidated method, the fragility function for a general IDR value is assumed to be a lognormal cumulative distribution function: where im is either S a (T 1 ) or PGA and Φ(·) is the standard normal cumulative distribution function. The values of the two parameters μ and σ, defining the fragility function, can be estimated using different techniques [60]. In the present paper the maximum likelihood estimation method was used. Considering the general IDR limit (either 1%, 3% or 5% in the present paper) and defining a vector Y containing the N values if im corresponding to this limit, the likelihood function for the observations in Y can be defined as: The estimates of the fragility model parameters (μ and σ) are those maximizing the likelihood function.
Fragility curves for the three IDR thresholds under consideration (representing slight damage, severe damage and collapse, respectively) were estimated both in terms of PGA and spectral acceleration S a (T 1 ). The median values of the fragility models, e µ , are reported in Tables 4  and 5 for the frames CM1 and CM2, respectively, and their logarithmic standard deviations (i.e. record to record variability) are reported in Tables 6 and 7. As an example, Fig. 8 shows some fragilities obtained for the structure CM2. In particular, two curves are plotted for each constitutive model: (i) the empirical cumulative distribution of structural capacity and (ii) the associated best fitting lognormal cumulative distribution. Fig. 8a shows results for IDR = 1% while Fig. 8b for IDR = 5%. Fig. 8, consistently with the results reported in Fig. 7, highlights the effects of the adoption of different constitutive models, both in terms of central value (i.e. median value of the fragility function) and in terms of variability. Considering the frame CM2 and IDR = 5% (Fig. 8b) the median values of S a (T 1 ), i.e. e µ , span between 0.336 g and 1.052 g, confirming the significant variability introduced by hysteretic models. This effect is relevant also for lower values of interstorey drift. For example, for frame CM2 and IDR = 1% (Fig. 8a) the median value of the curves ranges from 0.175 g to 0.259 g. It is worth noticing that the introduction of the in-cyclic strength degradation mechanism alone does not produce significant negative effects on the dynamic response of the structures. On the other hand, cyclic strength degradation mechanisms affect structural capacity. It should also be noticed that, as expected, record to record variability (σ) depends on the constitutive model adopted. Considering S a (T 1 ) as IM, the frame CM2 and IDR = 1% the values of σ range from 0.246 (PINCH-65) to 0.348 (IN-CYC_65), while for IDR = 5% from 0.305 (PINCH-65) to 0.502 (IN-CYC-65).
Then, it might be of interest to define global fragility curves for the structures, by considering all the structural capacity values obtained from the various constitutive models. These global empirical fragilities represent a more general set of curves. Fig. 9 shows the global empirical fragility curves obtained for the various structures in terms of spectral acceleration. Fig. 9a shows the results for IDR = 1% (slight damage condition) while Fig. 9b and c the fragilities for IDR = 3% (i.e. severe damage condition) and IDR = 5% (i.e. collapse condition), respectively. The median values (e µ ) and logarithmic standard deviation (σ) of the fitting lognormal distributions are summarized in Table 8, for both the intensity measures considered in the paper, i.e. PGA and Sa(T 1 ). The median values for the two structures is similar for IDR = 1%. While the logarithmic standard deviation ranges from 0.18 g (CM1) to 0.31 g (CM2). The next Section will show that for this IDR limit most of the standard deviation is due to record-to-record variability. Considering an IDR threshold of 3% the structures CM1 and CM2 show again a comparable median capacity. If compared with the data for IDR = 1%, total standard deviation is now significantly larger but similar for the two structures. Finally, for IDR = 5%, the structure CM2 designed by   Finally, as expected, it is worth noticing that the variability for fragilities in terms of PGA is higher than the dispersion of curves in terms of Sa(T 1 ), because of the higher efficiency of this latter ground motion intensity measure [59]. It is also interesting to highlight that increasing the IDR threshold, and therefore the extent of inelastic deformations, the total variance tends to increase. It should be noticed that, in the definition of the global fragility curves, every model has the same likelihood as being the true-model. As discussed in Section 6 it is possible to consider importance weights, which could be defined based on expert opinion, for the various models by associating weight coefficients to the different terms of the likelihood function in Eq. (4).

Analysis of variance
To evaluate the contributions to the total variance of the structural capacity values, related to record-to-record variability and model uncertainty, Analysis Of Variance (ANOVA) was used [61]. The ANOVA procedure used herein is based on a 2-way crossed classification model. According to this model the general observation of structural capacity y i,j , either in terms of PGA or Sa(T 1 ), can be written as: where μ is a general mean; i is the effect of the i-th ground-motion (i.e. y i,j is obtained using the ground-motion record i), j is the effect of the jth hysteretic model and i j , represents their interaction. The model in Eq. (5) is a random effects model which assumes that -,and -effects are random with zero means, variances 2 , 2 and 2 , respectively, and all covariance terms equal to zero. These assumption represent the customary formulation of random-effects in random or mixed models [61]. Assuming that the aforementioned factors are normal it is possible to write the distribution of the logarithmic structural capacity as: It should be noticed that Eq. (6) represents a simple analytical fragility model, similar to Eq. (3) [60,62]. The three components of the variance, for a balanced design of the numerical simulations (i.e. each ground-motion is used with each constitutive model), can be simply estimated as described in [56]. In particular, the analysis of variance table for the 2-way crossed classification model under consideration can be defined as shown in Table 9, where the different means are defined as follows: where g indicates the number of ground-motion records used, and c the number of constitutive relationships. The ANOVA estimators of variance components can then be computed as follows: The model parameters can also be estimated using the maximum likelihood method. The likelihood function, under normality assumptions, can be written as: where y indicates a vector containing the observed structural capacity values and V is the variance covariance matrix, defined as: The symbol indicates the Kronecker product, I and I are 11x11 and 30x30 identity matrices, respectively, and 1 , 1 11x1 and 30x1 unit vectors. The estimates of the parameters of the model, i.e. µ, , , 2 are the values that maximize the likelihood function (or its logarithm, which is easier to manage from a computational point of view). Because of the nature of the variance-covariance matrix which is sparse, specific computational methods are required for the solution of the maximization problem [61]. The model described so far can be further extended in order to consider importance weights for the various constitutive model. These weights could be defined based on expert opinion and engineering judgment. In order to include them in the fitting procedure, the easiest approach that can be followed is resampling structural capacity data proportionally to the importance weights, thus increasing the number of terms in the likelihood function. As an example, the weights reported in Table 10, based on engineering judgment, were used. It should be noticed that the judgement of different experts could be easily combined using logic tree.
Tables 11 and 12 report the results obtained for models CM1 and CM2, respectively. As expected, for low IDR values (1%) the model uncertainty has a marginal role because inelastic deformations are limited. On the contrary, for IDR values of 3% and 5%, the uncertainty related to constitutive model definition plays an important role, with variance contributions comparable to the record-to-record variability. For example considering IDR = 5% and Sa(T 1 ) these contributions represent 38% and 42% of the total variance for models CM1 and CM2, respectively. The contribution due to the interaction between model uncertainty and record-to-record variability is 8% for CM1 and 17% for CM2. This contribution represents the dependency of record to record variability on the constitutive model, as discussed in Section 5. For these levels of IDR, at which the structures show significant inelastic deformations, the choice of the constitutive model plays a crucial role in influencing the dynamic capacity of buildings. Therefore, when a   unique hysteretic constitutive model cannot be identified for a given structure, the introduction of the corresponding epistemic uncertainty, significantly increases the total variance of structural capacity.

Final remarks
The recent understanding of the importance of considering model uncertainties, especially in the definition of seismic collapse capacity has led the scientific community to focus their attention on this source of variability. In this context, the present paper mainly concerns the evaluation of the contribution to the total variability of fragility models, related to epistemic uncertainty on the choice of a proper hysteretic model. Two alternative designs, with and without capacity design rules, for a case study frame building, and eleven alternative constitutive models, based on the typical experimental behaviour of RC elements were considered in the study.
By means of IDAs we evaluated fragility for different IDR thresholds (i.e. 1%, 3% and 5%), considering either PGA or Sa(T 1 ) as groundmotion intensity measures. These data were analysed using ANOVA in order to identify and split the variance contribution due to record-torecord variability and model-to-model variability. We found that the variability contribution for model-to-model uncertainty is relevant for IDR values larger than 1% and it increases as the extent of inelastic deformation enlarges. At the collapse condition (IDR = 5%), magnitude of model-to-model variability and record-to-record variability are similar. Therefore, when the hysteretic constitutive model cannot be identified for a given structure, the introduction of the corresponding epistemic uncertainty, significantly increases the total variance of structural capacity.