Thick population density and its escalation propensity in seismically active regions of Pakistan has raised sincere concerns about the performance of building stock whose suboptimal performance and complete collapses led to a colossal number of casualties during the past earthquakes. The current research is inspired by the Kashmir earthquake of 2005 which consumed more than 80,000 lives, out of which, approximately 19,000 were children due to wide spread collapse of school buildings. A new database for existing reinforced concrete (RC) school buildings in seismic zone 4 of Pakistan has been developed using the surveyed information and presented briefly. The paper presents the statistics of the data collected through field surveys and professional interviews. It was found that the infrastructural authorities in the considered region developed some specific designs for school buildings, with varying architectural and structural configurations, which were eventually replicated throughout the area. In the current study, almost 2500 schools were surveyed for identifying versatile architectural and structural configurations, and subsequently, 19 different types had been identified, which were eventually used as representative stock for the schools in seismic zone 4 of Pakistan, Muzaffarabad district. The results of the study yield the brief of the collected data from the field and a consolidated methodology for establishing the analytical fragility relationships for one of the 19 structural configurations of the school buildings. A sample building from the collected data has been selected by considering the maximum number of students, and afterwards, the vulnerability is assessed by employing incremental dynamic analysis (IDA) which constitutes the presented methodology. Finally, the fragility curves are developed and presented for the said building type. The derived analytical fragility curves for the considered building type indicate its structural vulnerability and as a whole represent its satisfactory behavior. The vulnerability assessment process and the fragility development are described in an easy manner so that the domestic practicing engineers can readily become able to extend the application towards other school buildings in the region. The developed relationships can be employed for rational decision making so that essential disaster preparedness can be carried out by identifying any need for structural strengthening and interventions.

1. Introduction

The earthquake and its aftershocks that struck Pakistan in October 2005 were extremely catastrophic and served as an eye-opener for policy makers, engineers, and the hoi polloi on the necessity of having safe buildings and structures that must be resilient enough to sustain natural calamities. Seismically undermined infrastructural facilities and poor construction, specifically of school buildings, make the lives of youth severely vulnerable. Their safety and continuity of education is of paramount importance, not only for the parents but also for the long-term social growth and economic prosperity of the nation. A structure may get subjected to extreme seismic loading conditions during its life span, so an adequate engineering proficiency to assess diminishing structural capacity is essential to avoid complete structural collapse or to decide any specific intervention [1]. The 2005 Kashmir earthquake caused approximately 87000 casualties, out of which almost 19,000 were school going children, and cumulatively, it affected about 3.5 million people [2]. It was also stated by the Earthquake Engineering Research Institute Report [2] that approximately, 67% of the educational institutions in the area faced complete destruction. Despite of the deaths and socioeconomic ramifications, the 2005 earthquake yielded the phenomenon of “Build Back Better” in Pakistan, and the Building Code of Pakistan (BCP)–Seismic Provisions 2007 [3] is essentially a tangible example of it. In 2015, United Nations Development Programme (UNDP) reviewed the status of building design codes and bylaws with a precise focus on 2007 Building Code of Pakistan (BCP), and in its report [4], UNDP stated that a thin and tight development schedule forced a large-scale incorporation of American source documents and prevented the creation of a truly Pakistani code.

In major cities, developmental authorities with jurisdiction do most regulation over the respective cities and suburbs, but with respect to the BCP and earthquake design issues, the bylaws of leading local developmental authorities are incomplete and even inconsistent and at times, even contradictory [4]. From the prevailing information, it can be inferred that the design professional bears all the responsibility for structural design.

Considering the prevailing state of the vague practice of BCP 2007 and its provisions, the current research is targeted to assess the seismic vulnerability of RC school buildings that were designed and constructed in Pakistan’s seismic zone 4 after the 2005 Kashmir earthquake. Data collection was made in the whole of Muzaffarabad district of Kashmir as it is one of the most vulnerable zones and faced colossal devastation in the past due to earthquakes.

Numerous field visits were conducted, and professional interviews were taken with different organizations including National Engineering Services Pakistan (NESPAK), the most prominent organization to develop designs for most of the school buildings in the earthquake-prone regions. Acute effort was made to accurately develop the analytical model for vulnerability assessment. The field visits were conducted to visually inspect the contemporary condition of the school buildings and to take the essential measurements of structural members so that they may substantially be compared with the drawings obtained through professional interviews. Subsequently, using the incremental dynamic analysis (IDA), structural vulnerability is assessed and fragility curves are developed. As indicated by the literature, fragility curves are widely employed to represent the vulnerability of buildings [5, 6]. Tangible examples can be observed in the literature; for instance, Zain et al. [7] developed the analytical fragility curves for a reinforced concrete building in Philippines by considering two damage states. It is pertinent to mention that the damage itself characterizes undesirable changes in a system that adversely affects a structure’s behavior [8]. The development of fragility curves is not only limited to building structures. Pang et al. [9] and Rasheed et al. [10] conducted their researches for bridges, while Wang et al. [11] and Yazdani and Alembagheri [12] established the fragility curves for dams. The present paper specifically considers the building structure only for developing the analytical fragility curves. The analytical fragility curves can be generated using nonlinear static and dynamic analysis [13, 14], though nonlinear dynamic analysis is considered to be the most rigorous and reliable for their derivation [15, 16].

The present paper provides a brief insight into the collected data and presents a consolidated framework to assess the vulnerability of school buildings in seismic zone 4 of Pakistan. A newly developed database for the existing school buildings in the Muzaffarabad district of Azad Jammu and Kashmir (AJK) province of Pakistan is presented. For the purpose of developing analytical fragility relationships, a building typology, representative of the typical school building stock with certain structural configuration, is selected for developing its fragility relationships against seismic loading. Though some work has been done in this sphere in the developed world, but developing countries still lack such researches, and presently, no other work exists that can specifically address the development of analytical fragility relationships for school buildings in Pakistan. The outcomes of the presented study can be readily adopted by relevant authorities and disaster preparedness and mitigation agencies to decide any need for specific structural or nonstructural intervention. The application of the presented framework can be further extended to establish analytical fragility relationships for other types of school buildings whose information are briefly presented in this work.

2. Data Collection and Methodology

2.1. Data Collection

Data collection was made from whole of Muzaffarabad district to establish a database for existing reinforced concrete (RC) school buildings. The database contains all the information pertaining to structural features of schools, i.e., number of stories, locations of staircases, locations and dimensions of beams and columns, thickness of slabs, and the width and locations of infills. The database also contains the structural drawings, i.e., design and as built that were collected as part of the data collection process. Later on, the obtained drawings were compared with the existing structures to obtain insight into the actual execution and the contemporary condition of the school building structures.

2.2. Development of Database for Existing School Buildings in Muzaffarabad District

A total of 2417 schools were visited, and visual inspection was made in accordance with FEMA 154 [17]. Detailed plans, indicating the locations and sizes of structural components, were developed for all the buildings and were eventually compared with the obtained drawings from different domestic sources. Out of 2417 schools, entrance was granted for 1933 school buildings, and the local staff for measuring the dimensions of the structural components extended appropriate cooperation and elaborated about the behavior of the buildings during the earthquakes which came after 2005. While for the rest of the buildings, i.e., 484, only exterior was available for observance from which the number of bays, external columns, and beams could be observed. It was found that out of 2417 school buildings, 1158 were single story, 847 were double story, and 412 were three-story structures. It can be conveniently deciphered from the collected data that most of the building population in the region consists of single story school buildings, with obvious variation in the number of their bays and story heights.

It was found out that domestic authorities developed some typical architectural and structural designs which were subsequently replicated throughout the considered region for constructing RC schools. For the purpose of developing categories of school buildings with varying structural configurations to assess their vulnerability, 19 different types of RC school buildings were identified through the field surveys and professional interviews. The inventory developed after the identification of the building types, along with their structural features, is given in Table 1. The first column of the Table 1 denotes the name of the model, aimed to represent the generic nonlinear model for that category of schools. All the names start from “BLR,” and subsequently, the model number is mentioned, i.e., BLR-1, BLR-2, up to BLR-19, representing all the 19 school building typologies with varying structural and architectural configurations. “BLR” in the names means “Building-Low-Rise,” essentially representing the low-rise nature of the construction in the considered region as the number of stories for school buildings remains between 1 and 3 only. The 2nd column presents the number of stories; column 3 and 4 provide the number of bays in x and y directions, respectively. The x and y directions are referred as the direction along the front side and the orthogonal direction, respectively. It was observed that mostly two bays existed in the y-direction, i.e., one for the pedestrian walkway which mainly served as a corridor so that students could walk to get in and out of the classrooms, while the other one served the purpose of a classroom. The total height of the buildings is given in column 5. Column 6 presents the typical floor area while column 7 corresponds with the area of slab openings, placed above the staircases. Typical floor area means the area for a single floor only.

Figure 1 provides the number of schools under each of the 19 school building categories. It was found that BLR-5 was the mostly constructed school building in one-story structures; BLR-11 was the mostly constructed school building in two-story structures; while BLR-18 was the mostly constructed school building in three-story structures in the considered region.

In the present study, the vulnerability is assessed by considering, a two-story structure, BLR-11. The BLR-11 is specifically considered as this typology has the maximum number of students as per the collected information. The building categories with 3 number of stories housed allied facilities, i.e., libraries and laboratories; but it was found that among all typologies, BLR-11 accommodates most of the pupils from the considered seismic zone. The altruistic intention is to produce building class-specific fragility relationships that can be easily adopted by the researchers and, more practically, by practicing engineers. The fragility curves have been developed using the IDA, the details of which have been provided under the subsequent section.

2.3. Methodology

For conducting the vulnerability assessment of school buildings in seismic zone 4 of Pakistan, the current study employs nonlinear incremental dynamic analysis (IDA) to capture the structural response against applied and incrementally scaled set of ground motions. For this purpose, 20 ground motions were carefully selected and scaled at increments of 0.20 g, starting from 0.20 g up to 1.40 g.

Since earthquake loading parameters are usually characterized by single numerical values, i.e., peak ground acceleration (PGA) or pseudo spectra acceleration at the fundamental period of the building (T1), single-value indicators, PGA and Sa, at the fundamental period of the building (T1), have also been used in the current study to correlate with the structural damage. As low-rise buildings are the primary focus of this research, the current study employs global drift to decide the threshold damage limits to decide different damage states, also known as limit states, of the building. Dynamic instability is usually caused by large global drift due to large interstory drift and P-delta effect. Keeping this in view, global drift is selected as the damage indicator to incorporate the overall structural global stability in the presented work.

Figure 2 qualitatively depicts the illustration of the deformational-based response, corresponding with different limit states used by the HAZUS [18]. The figure illustrates four damage states, ranging from slight damage to complete collapse of the building, by relating the base shear with the deformational-based structural response, while the deformational-based response itself can be considered in numerous forms, i.e., global drift, interstory drift, or the spectral displacements.

Finally, using a probabilistic seismic demand model attained by the virtue of simulated damage data, fragility curves for each of the considered damage state were obtained. The framework and procedure through which the vulnerability assessment has been conducted is presented in Figure 3. Detail for every step is given hereafter.

3. Application of the Framework

For demonstrating the established framework, a nonlinear model for a selected school building type, BLR-11, from the developed database is created using CSI Perform 3D V 7.0 [19] as it contains the vital inelastic elements that can monitor the strains in structural members. For the considered typology, nonlinear inelastic fiber elements were used to monitor the strains in the members and substantially, to set threshold limits for the attainment of damage states. The configuration of the selected building, employed analytical tools, and material strengths are discussed in the forthcoming sections.

3.1. Configuration of the Selected Building Typology

The selected school building represents a typical two-story structure that comprises 11 bays in the x-direction and 2 bays in the y-direction. It consists of different sizes of beams and columns. Figure 4 shows the plan of the selected building that has been considered in the present study, while Figure 5 exhibits the elevation of the considered building in the longer dimension to indicate the story height.

In Figure 4, FB means the floor beams, while FC means the floor columns. The actual buildings are constructed on relatively stiff soil type, but it was found during the professional interviews that all of them have been designed according to the soft soil type, Sd soil according to BCP. The main structural system of the buildings consists of moment-resisting RC frames in both of the directions. The typical floor area for the selected building is 344.5 m2. Table 2 shows the X-sectional dimensions of structural members, i.e., different FBs and FCs. The amount of reinforcement, obtained through field surveys and professional interviews, is also indicated in the same table.

3.2. Nonlinear Structural Modelling

Full three-dimensional nonlinear structural model was created in software CSI Perform 3D v 7.0. The analytical tools in the software possess excellent capability to capture different aspects of structural performance, ranging from hinge rotations to materials’ strains in members. In order to capture the actual behavior, it was pertinent to take into account the properties of locally developed concrete materials using the domestically available cement, aggregates, and other constituents. For this purpose, a study conducted by Rafi and Nasir [20] was utilized to extract a mean value for concrete’s strength at 28 days, subsequently, as per the recommendations of ASCE 41 [21], expected strength was ultimately employed during the analysis. For reinforcing steel bars, it was inferred from most of the professional interviews that Grade 40 reinforcing bars were used during the design and construction process, consequently, the current study uses Grade 40 steel reinforcements in the analytical models in accordance with the recommendations of ASCE 41 [21].

Nonlinear inelastic fiber sections were created for layer-by-layer modelling of concrete and reinforcing rebars for all of the X-sections except slabs, which were modelled as linear. By the virtue of fiber modelling, the force-deformation relationship of columns and beams is transformed into the stress-strain relationship of construction materials, i.e., concrete and reinforcement. By employing fibers, an appropriate division may be introduced for confined and unconfined concrete fibers and consequently, a composite section can easily be taken into account [22].

In the current study, the Mander model, provided by Mander et al. [23], was used for concrete. Complete confinement effect and strength loss were considered during the analysis according to the details given in the X-sections of beams and columns (Table 2).

Steel models that can capture the buckling and nonbuckling behavior of the reinforcements are available within the Perform 3D. For this case, the nonbuckling steel model was utilized as ductility design is mainly based on the fact that reinforcement cannot be abruptly brittle. Expected material strength was used in the analysis as per the recommendations of ASCE 41 [21]. The numerical modelling of strength degradation at high ductility levels was also considered by modelling the strengths of the materials in a nonlinear range. Table 3 summarizes all the material strengths, live loads, and infill-partitioning loads that were used during the analysis in accordance with the professional interviews conducted during the field surveys.

4. Uncertainty Treatment

The uncertainties can be divided into two categories, i.e., aleatory and epistemic during any sort of vulnerability and risk assessments [24]. The uncertainty associated with the inherent variability in the nature of ground motion, resulting from the natural earthquakes, comes under the heading of aleatory uncertainty. Thus, the intrinsic randomness and variability of ground motion record describe such uncertainties.

The epistemic uncertainties, on the other hand, characterize the qualms related to paucity of available knowledge. Conventionally, during the vulnerability assessments, the lack of knowledge related to the materials’ strengths comes under the category of such uncertainty. On the contrary, Zain et al. [25] have shown that variation in material characteristics do not induce any significant disparity in the structural response, and variation in the ground motions must be considered as the majorly significant factor in influencing the dynamic response of the structures, as record-to-record variability may persuade severe changes in the structural response of any building against the properties of different ground motions.

The current study employs 20 different ground motions to cover a wide range of variation in the characteristics of seismic loading. Since the variations in the materials’ strengths do not incite any major uncertain structural response, hence their variation is not considered in this study. The details about the selection procedure and incremental scaling for the ground motions are discussed under the proceeding section.

4.1. Ground Motion Uncertainty

The biggest source of uncertainty lies in the seismic demand during the vulnerability assessment process. Path attenuation, characteristics of the originating source, and the site conditions all play their respective parts in contributing to such uncertainty. On account of the structural dynamic response, it is imperative to consider a wide range of seismic energy levels for selecting the ground motions. Typically, a target spectrum is considered for the anticipated seismic hazard level, and the spectrum-matching technique is usually adopted to select the most suitable ground motions. Although, the Building Code of Pakistan (referred at [3]) provides a target spectrum for the region considered, seismic zone 4 of Pakistan; however, as explained in Section 1, the reliability of the domestic code is questionable. Therefore, in the present study, ground motions are selected by considering the best available information of the case study area. In this research, 20 natural ground motion records are selected by considering the domestic fault mechanism, magnitude, source-to-site distance, and the shear wave velocity Vs, in the upper 100 feet of the soil layers, and only the large magnitude earthquakes, i.e., 6.50 up until 8.0, are considered.

For the current research, the range of source-to-site distance has been kept between 0 and 30 kilometers for selecting the ground motions to consider the effects of path attenuation so that any possibility of diversified site-specific soil layers in near and far sites can be taken into account. For shear wave velocity Vs in the upper 30 meters of the soil, the range has been kept varying between 175 and 350 m/sec which corresponds with the soil type Sd according to BCP-2007 (referred at [3]).

The selected ground motions are given in the Table 4, which presents the names of earthquake motion time histories, magnitudes, distances to rupture, shear wave velocities, and associated PGA with some other minor details.

By using numerous records, with the consideration of domestic geological conditions, the selected records may be considered as the representative of future shaking at the zone in question. Thus, with a variety of magnitudes, PGA, and other characteristics, the aleatory uncertainty is adequately treated in the present study.

5. Incremental Dynamic Analysis and Ground Excitation Random Effects

Incremental dynamic analysis (IDA) serves as the most reliable method to assess the structural performance against seismic excitations [26]. Several researchers, i.e., Kostinakis and Athanatopoulou [27], Fereshtehnejad et al. [28], and Zarfam and Mofid [29] adopted the technique for the evaluation of dynamic responses of different structures, i.e., buildings and bridges, and despite of its extensive computational effort, IDA has proved to be a splendid tool for probabilistic vulnerability and risk assessments. During IDA, the ground motions are incrementally scaled in a way that during each iteration of the analysis, the accelerations increase at specific intended intervals. The current study utilizes the IDA for considering various levels of seismic intensity to predict the structural response against earthquakes, ranging from 0.2 g to 1.4 g with 0.20 g as the augmenting acceleration for each successive iteration, thus, at this stage, it is imperative to select the appropriate seismic intensity measure (IM) for representation of seismic. Subsequent subsection addresses the selection of IM in the current study. Apart from selecting the IM, there must also be a damage indicator that can amply relate the seismic IM with the observed structural damage. As described earlier, the current study employs global drift as the damage indicator, and the details have been given in Section 6.1.

5.1. Seismic Intensity Measure

Seismic intensity measure (IM) correlates the ground motion with the structural damage; hence, any IM employed must be rationally able to project the dynamic response of the structure by the virtue of the damage indicator. PGA is one of the IMs that has been frequently used by the experts for fragility derivations, but it is generally comprehended that Sa is a better option as an IM to correlate with the structural damage as it actually represents the effect of ground shaking on the building itself. Pang and Wu [30], and Jiang et al. [31] employed PGA for their researches. Apart from the PGA, Omine et al. [32] and Eden et al. [33] used peak ground velocity (PGV) as an IM. Eden et al. [33] also used Sa at the fundamental period of the structure as an IM. Frankie et al. [34] and Choudhury and Kaushik [35] utilized spectral displacement (Sd) as IM for their study on open ground story RC frames with openings in walls during the vulnerability assessment. Buratti [36] analyzed the sufficiency and the efficiency of numerous IMs as indicated in their study and employed the cloud analysis technique for 1704 unscaled accelerograms. Bojorquez et al. [37] proposed a procedure to obtain the interstory drift demands for midrise frame structures in terms of a spectral-shape-based IM () and eventually concluded that was the best parameter to relate with the structural response of frames with 4, 6, 8, and 10 stories under the narrow-band ground motions in Europe. Some other vector-valued intensity measures have also been proposed and discussed by Modica and Stafford [38], Baker and Cornell [39], Polsak and Nicolas [40], and Kohrangi et al. [41]. Bojorquez and Iervolino [42] also performed their research to establish a vector-based IM; however, their research also proposed a scalar-based IM as they apprehended that the use of the scalar-based IM is easier than the vector-based IM for practical purposes.

Pejovic and Jankovic [43] investigated different ground motion IMs which are conventionally used to investigate the dynamic behavior and concluded that despite being the most common IM, PGA provided the greatest dispersion in results in comparison with the PGV, Sa(T1), Sv(T1), and Sd (T1). The research portrayed spectral response velocity value Sv(T1) as the most efficient IM for the RC frames, while stating the other two spectral response parameters, i.e., Sa(T1) and Sd(T1) were similarly adequate as they incorporated the dynamic characteristics of the considered structure. The current study chooses both, i.e., PGA and spectral acceleration at the fundamental time period (Sa(T1)) of the structure as intensity measures. As indicated by the previous literature, provided by Jelena and Srdjan [43], PGA is the most commonly employed IM which is easy to understand for normal people with no specific engineering background, while Sa(T1) incorporates the structural characteristics during the execution of analyses.

5.2. Results of IDA

Figure 6 shows the results of IDA in the form of plots between global drift and the selected IMs, i.e., incrementally scaled PGA and Sa(T1) from 0.20 g to 1.4 g. It can be observed that intrinsic nature of different ground motions produces different results despite of having the same intensity of acceleration at a given point. There are total 20 curves in each part of Figure 6, while each curve characterizes structural damage, indicated by global drift, caused by the single corresponding earthquake at discretely considered seismic intensity levels, ranging from 0.2 g to 1.4 g.

In Figure 6(a), the curves manifest the maximum global drift ranging from 0.20% to more than 5.0% when PGA is used to indicate seismic intensity. While in Figure 6(b), when Sa(T1) is used to characterize the seismic intensity, the maximum global drift reaches 2.47% against the Sa of 1.4 g, the maximum considered in this study.

The results of IDA indicate that intrinsic properties of ground motions project substantial sensitivities to the seismic demand. The disparity in the structural response against different ground motions of the same intensities is primarily attributed to differences in the inherent nature of ground motions, i.e., energy distribution over frequency content and duration that predominantly affects the structural response. Moreover, when Sa(T1) is used as an indicator of the seismic intensity, the global drift observes lower values because of the scaling of the ground motions at the fundamental period of the structure which is directly related with the structural inertia and stiffness.

6. Vulnerability Assessment

The fundamental process of the risk assessment against any hazard consists of two parts, i.e., the hazard and the vulnerability of structural systems to that hazard. For earthquakes, the first portion is primarily concerned with seismologists who can provide a realistic estimate of the potential seismic hazard that can hit a specific geographical area; while the second part, vulnerability, is usually addressed by the engineers and experts who design the facilities for versatile purposes. The current study is related with the second part and assesses the vulnerability of the considered typology of school buildings by establishing the hazard-damage relationships for the considered typology of structures and subsequently, by developing the generic fragility relationships. Hazard-damage relationships depict the structural performance against seismic excitations and exhibit the dynamic response for all considered ground motion intensities. Figures 7(a) and 7(b) present a direct relationship between the seismic intensities and the structural damage. The structural damage is quantified using the global response, and seismic intensities are provided incrementally in terms of PGA, as well as in terms of spectral acceleration, at the fundamental mode’s time period. Every single dot in each part of Figure 7 represents the maximum value of the global drift, in terms of percentage, against a specific earthquake, scaled at a specific intensity. For instance, at the seismic intensity of 1.4 g, there are 20 dots from all 20 ground motions, indicating the uncertainty in the structural response by depicting that one earthquake having the seismic intensity of 1.4 g produces global drift less than 1%, while some other earthquake with the same seismic intensity produces more than 5% of global drift.

With increasing intensity of PGA, from 0.20 g to 0.80 g, in Figure 7(a), the global drift remains between 0.09% and 1.19%. While at higher PGA levels, from 1.0 g to 1.4 g, the global drift drastically remains between 0.44% and 5.16% for all the considered ground motions. At the maximum considered PGA of 1.40 g, the global drift varies from 0.97% to 5.16% for all 20 ground motions. All this disparity in the structural response arises because of the inherent characteristics of ground motions which vary from record to record, thus instigating different response against the same scaled intensity level.

When Sa(T1) is used for representing seismic intensity, the global drift varied between 0.22% and 2.4% at the Sa of 1.40 g, as shown by Figure 7(b). The dotted line in Figures 7(a) and 7(b) shows the mean structural response obtained at each intensity level.

On the other hand, the fragility curve, as described by Equation (1), represents a conditional probability of exceedance for specific damage or limit states of a structure against a particular seismic intensity given in terms of PGA, Sa, or some other seismic intensity describing the indicator.

The aforementioned equation describes the conditional probability, Pfragility, of achieving or exceeding a specific limit state, LS, when seismic intensity, defined by a particular intensity measure (IM), is equal to “x”. Fragility relationships for different damage states must be probabilistically able to predict the structural vulnerability against a realistic range of seismic ground excitations. This probabilistic evaluation is only possible when a large variation of seismic demands is considered so that ample variation in the dynamic structural response can be obtained. The vulnerability assessment requires the definition of specific limit states, damage indicators, and seismic intensity indicators. The impending sections address all the requirements to develop fragility curves. Consequently, fragility curves are developed for the selected typology to elucidate the established framework.

6.1. Damage Indicator and the Definition of Limit States

The derivation of analytical fragility curves requires the definition of specific limit states in both terms, i.e., qualitative and quantitative. For this purpose, a researcher has to employ a specific damage indicator, also known as engineering demand parameter (EDP), or may also propose a new damage measure that can represent the structural damage in a rational manner to correlate the damage with seismic loading. Once the damage indicator has been set or selected, different limit states, often known as the damage states, of a structure may be defined.

Damage states represent a physical categorization of damages by using damage indices which suggest threshold limits for different damage states by relaying upon the parameters of structural behavior, e.g., energy dissipation capability, ductility, and strength. Damage indices are used for correlating the accumulated damage with different specific damage states, and numerous nonidentical damage indices may exist for the same structure. For instance, Khan et al. [44] used global displacement of the structure to characterize the structural damage. Crowley et al. [45] considered damage limit states based on the strain levels in concrete and steel. Güneyisi and Altay [46] considered the interstory drift ratio for representing the response parameter corresponding to the four damage states, taken from the HAZUS. Casotto et al. [47] considered the member flexural strength to characterize the first limit state (LS1) when the reinforcement steel yields in the columns and 3% interstory drift for flexural collapse limit state (LS2), and subsequently, related his LS2 with the loss of support of the beam. Chaulagain et al. [48] and Gautam et al. [49] used 3 damage states ranging from slight damage to complete collapse and developed empirical fragility relationships for buildings in Nepal. They further related their results with FEMA356 [50] classifications for building damage levels.

Since, this study primarily relates with the generation of fragilities for a specific building stock, it would not be practical to define damage or limit states depending upon the member behavior or local strains. In the current study, the global drift ratio has been considered as the damage indicator to correlate with the structural damage.

Three damage states, namely, serviceability limit state (LS1), damage control limit state (LS2), and collapse prevention (LS3), are qualitatively defined in this study. LS1 corresponds with the 1st yield of longitudinal steel reinforcement in the structure. Conventionally, the first limit states are usually related with the formation of the first hinge in the structure, but since the current study employs nonlinear fibers to model the structure, therefore, the LS1 primarily describes the yielding strain reaching in a structural element, located anywhere in the structure. The deformation and strength govern the LS2. It is defined as 75% of LS3. The LS3, collapse prevention, is generally directed by the deformation. Altug [51] defined the collapse prevention limit state as 75% of ultimate structural deformation or the value for which more than 20% strength drop can be observed relative to the maximum strength value. The smaller of these two values would govern the LS3.

The mathematical quantification for these damage states could not be adopted from the previous research as there exists no analogous research for Pakistan’s construction industry or for Pakistan’s seismic zone 4. In the current study, the mathematical values are assigned to each of the discrete damage state by using the capacity curve of the selected structure. However, the obtained values were approximately similar to those obtained by Altug [51], which used the low-rise buildings of the Turkish region in their research, thus the obtained values amply represent the low-rise structural behavior. The mathematical values for each of the considered damage state against the global drift are presented in Table 5.

6.2. Fragility Assessment for the School Buildings

Seismic fragility curves are widely employed to obtain an articulated insight of structural performance against seismic loads, and they represent the probability of entering or exceeding different damage states for considered intensities of ground shaking. After attaining the structural response, established limit state definitions are applied for evaluating the structural performance for each limit state. Once the computed value comes out to be greater than a limit state’s value, an event shall be counted within the sample, comprising all the ground motions scales. Thereby, for each IM, direct sampling probabilities over the 20 selected ground motions on all incrementing scales were obtained.

It is a commonly accepted assumption that for a fragility curve, a lognormal cumulative distribution function can be conveniently utilized for the representation of conditional probabilities [5254]. From viewpoint of literature, lognormal distribution was also assumed and used for the regression of fragility relationships as mentioned in the following equation.where describes the probability of exceedance of a particular LS, given an IM value. represents the standard normal cumulative distribution function (CDF).

The controlling parameters and represent the median point and the slope of the curve, respectively. It is imperative to evaluate their optimized values for a best fit of CDF against the calculated probabilities. Shinozuka et al. [53], Dang et al. [54], and Baker [55] stated that the maximum likelihood method (MLM) is one of the most suitable method to achieve best possible values of these two parameters. In the current study, MLM is utilized to obtain the median point and the slope of the fragility curve. For multiple IM levels, as considered in this study, the likelihood function is given as follows:where m is the number of IM levels and Π denotes a product over all levels. It is assumed that the observation of attaining a specific limit state from each ground motion is independent of the observations from other ground motions. In the above equation, zj denotes the number of attaining particular limit state out of nj ground motions with a particular value of IM = IMi.The fragility curve parameters can be obtained by maximizing the likelihood function in Equation (3) as given in following equation.

A MATLAB code was written to control the maximization of the likelihood function so that optimum values of these two parameters can be obtained for establishing seismic fragility.

Figure 8 presents the developed fragility relationships for the considered typology of the RC school building in seismic zone 4 of Pakistan. Each portion of Figure 8 contains three fragility relationships correlating with each limit state against both of the IMs, i.e., PGA, and Sa(T1). Direct sampling probabilities, obtained from the IDA, are also plotted in Figure 8 along with lognormally fitted fragility functions.

Table 6 provides the simulated values for and for each limit state against each IM for the derivation of fragility curves.

7. Conclusion

Pakistan has numerous seismotectonically active regions including Muzaffarabad and past earthquakes, especially the 2005 Kashmir earthquake, have proven this fact. In this study, a new database of RC school buildings, constructed in Muzaffarabad district after 2005 earthquake, has been presented. Subsequently, the structural vulnerability of a specific category of RC school buildings is assessed and presented by using IDA so that practicing engineers at domestic and regional level can further extend the evaluation to other categories of RC buildings. The procedure presented utilizes the fiber-based model which is the most reliable at present. In the current study, both PGA and Sa(T1) have been used to characterize the seismic intensity. Though Sa(T1) is a more appropriate IM, PGA has been easily understandable with the normal public; therefore, the fragility relationships have also been developed using the PGA as an IM. At the PGA of 0.6 g, there exists approximately 82% probability of exceedance for the first limit state (LS1), while only around 20% probability of exceedance has been achieved against the same seismic intensity for the second limit state (LS2). This essentially indicates the structural behavior dominated by the yielding of the members instead of their complete failures. At the seismic intensity of 1.40 g of PGA, the structure seems to be in LS3 with almost around 100% of probability; while when Sa(T1) is used as an IM, the probability of being in LS3 at the seismic intensity of 1.40 g remains only around 20%. This specific variation in structural behavior is primarily attributed to the inherent structural characteristics, i.e., stiffness and ductility that influence the response when the spectral acceleration is used as an IM, after being scaled at the fundamental time period of the structure. Since, the selected structure is engineered, the results are not so admonishing for the considered school building typology; nevertheless, the results might not be encouraging for other school building types in the presented database. The presented consolidated procedure is relatively simple and can be readily adopted by the domestic and regional disaster preparedness and mitigation agencies to extend the application. Such studies are relatively new for under developing countries like Pakistan, and contemporarily, no other work exists that may specifically address the vulnerability of schools in seismic zone 4 of Pakistan. Assessment of vulnerabilities for other types of school buildings may yield a need for immediate intervention in order to safeguard the life of students and other users. The presented study is specifically targeted to assess the vulnerability and does not address the assessment of risk; therefore, the developed fragility curves can be further integrated with the hazard curves of the considered area to evaluate the total risk involved.

Data Availability

The data used to support the findings of this study are included within the article. Any other data used to support the findings of this study may be released upon request to the author(s), who can be contacted at [email protected].

Conflicts of Interest

The authors declare that they have no conflicts of interest.