Abstract

The paper presents the development of a nonlinear static displacement-based methodology for seismic risk assessment and loss estimation of stone masonry building stock of Pakistan. Experimental investigation of one-third scaled model, tested on shake table, is performed in order to obtain lateral strength and drift limits for stone masonry and develop damage scale for performance-based assessment. Prototype buildings are designed respecting the existing building stock and investigated through nonlinear static and dynamic time history analyses. Nonlinear static mechanical models, for both global and local vulnerabilities, are developed for the considered typology which are used to derive analytical structure-dependent fragility functions considering expected sources of uncertainties explicitly in contrary to the conventional procedures. Furthermore, seismic risk assessment is performed for different scenario earthquakes and presented in terms of structure-independent fragility functions to estimate the mean damage ratio, the repair cost as a fraction of replacement cost, and casualties, with the dispersion being quantified, given source-to-site distance and magnitude for an earthquake event. The methodology is tested for seismic risk assessment of the considered typology in recent 2005 Kashmir earthquake, which is reasonably predicted. Future development of the methodology is required with additional experimental tests on rubble stone masonry material in order to increase confidence in future applications.

1. Introduction

Stone masonry buildings constitute a substantial portion of urban and rural building stock of northern areas of Pakistan. Two wythes random rubble stone masonry walls (Figure 1) in dry or mud mortar or cement mortar with flat earthen roof (wooden beams and straw), wooden/steel truss with galvanized iron (GI) sheet roof, and reinforced concrete (rc) slab are the most common construction type in rural and urban settlements of Pakistan. Recent practices in the country, mainly in the recent earthquake affected areas (http://www.erra.pk/), make use of throw stones and vertical steel bars, 1.2 m apart, in order to improve the seismic performance of such systems which seems not to be dramatically different than the ordinary wall configuration. Nevertheless, the findings herein are conservative for the new structural systems. Furthermore, the overall contribution of old construction schemes is relatively more in Pakistan [1].

Such building systems have shown very poor performance in past earthquakes and lead to huge losses of life and economy. Collapse of such structures featured prominently in the 2005 Kashmir earthquake [2, 3]. However, due to the abundant local availability of stone material, cheap labor cost, low economic status of the people, old traditions of construction, and lack of public awareness, stone masonry buildings are still in use in Pakistan. Stone masonry buildings are practiced, with little modifications (Figure 1), in the earthquake affected areas in the order of 50 percent of total building stock, which is about 500,000 housing units. Also, some parts of the country, which are not affected in the near past, have significant amount of old stone masonry buildings (Figure 1) [1]. Thus, the present study assesses the seismic performance of the considered building system using experimental and numerical investigations. Analytical methodology is developed for future applications in seismic risk assessment and loss estimation of stone masonry buildings on regional scale. Simplified charts are developed to estimate the socioeconomic impacts of earthquakes with the expected dispersion being quantified, given the magnitude and source-to-site distance of event, for earthquake preparedness and planning, decision-making, and risk mitigation in Pakistan. A test prediction of seismic risk is performed for the considered typology in the recent 2005 Kashmir earthquake which is reasonably predicted.

2. Displacement-Based Earthquake Loss Assessment Methodology (DBELA)

Seismic risk assessment and loss estimation of structures have developed from initial studies, in the past 30 years, which used empirical approaches to predict the socioeconomic impacts of earthquakes, to more recent methodologies based on numerical predictions of the capacity of buildings to resist strong ground shaking [4, 5]. There are various shortcomings with the use of empirical methodologies, which arise due to the use of ground motion intensity parameters, which do not take the relationship between the frequency content of the ground motion and the vibration characteristics of the building stock into account (e.g., peak ground acceleration). The ground motion uncertainties and spatial variabilities cannot be truly represented which can be different from earthquake-to-earthquake and site-to-site. Further problems exist with such methods due to the use of existing building damage, which thus precludes the possibility of accurately modeling the vulnerability of new or retrofitted buildings. In order to respond to these shortcomings and, perhaps more importantly, so as to allow the most recent developments in seismic hazard and ground motion prediction to be employed, newly emerging methodologies are based on the analytical assessment of building damage given an input response spectrum. Such methods have the additional advantage of allowing the vulnerability of the building stock to be predicted based on a limited amount of knowledge of the structural characteristics of the buildings under consideration and thus they can be easily calibrated around the world [5]. Furthermore, the vulnerability predictions with mechanics-based methods can be easily justified following a close inspection of the characteristics of the buildings. Thus, the present study performed calibration of a nonlinear static analytical mechanics-based, mainly displacement-based, fully probabilistic method, presented herein, for seismic vulnerability and risk assessment of stone masonry building stock of Pakistan on regional scale. The displacement-based method is originally proposed and developed elsewhere for rc and masonry buildings [610]. However, it is further extended and developed to take different sources of uncertainties, involved in the evaluation of structural capacity and seismic demand, in an explicit and transparent manner and make use of analytical fragility functions, derived in a conceptual and transparent manner, to develop damage scenario for a given earthquake. Also, the methodology can be applied to both individual building and building aggregates considering their global and local vulnerabilities. The methodology is previously developed for other building stock of Pakistan [1115] and further calibrated for stone masonry buildings of Pakistan herein.

2.1. Fundamentals of the Methodology

The methodology makes use of analytical fragility functions and a mathematical model for building stock, which is represented as a nonlinear static single degree of freedom (SDOF) system defined completely by secant vibration period, displacement capacity, and viscous damping, to assess the seismic risk and socioeconomic impacts of a given earthquake event. The seismic demand on building stock is defined by 5% damped displacement response spectrum using code-based (for intensity-based), GMPE-based (for scenario earthquakes), and/or UHS-based (for estimation of loss exceedance curve and annualized losses) spectrum representation depending on the scope of the study. GMPE stands for ground motion prediction equation and UHS stands for uniform hazard spectra. To better understand the methodology, it is depicted graphically (Figure 2).

In this methodology, the first step is to generate elastic displacement response spectrum for a given earthquake, which is overdamped in an iterative fashion to compute the structure-specific seismic demand (seismic intensity). This seismic intensity is used to compute the number of buildings in different damage states from the analytical fragility functions in order to develop damage scenarios. The number of buildings in different damage states is transformed to mean damage ratio (MDR), which quantifies the regional seismic risk and gives an estimate of direct economic losses for the considered earthquake. The conversion of building damage states to MDR is performed using earthquake loss model [16]. However such models are needed to be calibrated for different regions in order to develop regional specific risk estimation tools and increase the confidence in future applications.

2.2.  Mechanical Model for Stone Masonry Buildings

For seismic assessment, the stone masonry building stock of Pakistan is broadly classified in two major classes based on the predominant seismic response mechanism during an earthquake. Masonry buildings provide seismic resistance to ground shaking through the development of in-plane/out-of-plane forces/deformations in structural walls. The buildings that have enough structural integrity provide seismic resistance predominantly through in-plane shear/flexure response of structural walls (global mechanism) while buildings with less structural integrity respond primarily in out-of-plane response of walls (local mechanism). The methodology uses SDOF systems to assess the seismic performance of each class of building stock. The SDOF system, called a mechanical model, has nonlinear lateral force-displacement response to assess the seismic performance of structures. The mechanical model simulates the response of the structural system in terms of its displacement capacity, energy dissipation, and secant vibration period for seismic assessment. The SDOF system derivation for each class of buildings are discussed as follows.

2.2.1. Global Mechanism (In-Plane Response)

The seismic response of stone masonry buildings having rc slab, well-connected orthogonal walls, and ring beams is mainly governed by the global response of buildings, because of the in-plane integrity of walls provided by the rc slab, with shear failure dominated mechanism of in-plane walls [1719]. Flexure rocking with compressed toe crushing is also a possible mechanism to resist lateral load. The structural system provides energy dissipation due to nonlinear behavior of masonry walls and limited ductility capacity. The current construction practice in the area is to include vertical steel bars and through stones in wall; although it may not increase the wall strength significantly, however, it will ensure sufficient in-plane integrity of the walls and whole system to respond through global mechanism. An equivalent SDOF system (Figure 3) is used to simulate the global response of building in terms of displacement capacity, viscous damping, and secant vibration period at different damage states for performance evaluation given an earthquake event, represented as 5% damped displacement response spectrum.

In this figure, represents the total building height; represents the th floor height, represents the lateral displacement, and represents the th floor mass for a given deformed shape of building; and represent mass and height of the equivalent SDOF system; and represent the equivalent yield and ultimate limit states displacements that represent the displacement capacity of the actual building at the center of seismic force for a specified deformed shape; represents the initial preyield stiffness; represents the yielding force; represents the secant stiffness.

For seismic assessment, the mechanical model is completely defined by secant vibration period, limit state displacement capacity, and energy dissipation characteristics of building represented as viscous damping:where represents the limit state secant vibration period; represents the yield vibration period; μ = represents the limit state ductility; represents the yield displacement capacity; represents the specified limit state displacement capacity; represents the yield drift; represents the specified limit state drift; represents the height of the system; and represent the displacement coefficients to convert multi-degree of freedom (MDOF) structural system to an equivalent SDOF system and simulate the displacement capacity of MDOF system at the center of seismic force [8]; represents the equivalent viscous damping of the structural system [20]; represents the elastic damping of the system (preyield); represents the hysteretic contribution of system damping due to nonlinear response of the structural components. The limit state parameter’s values are selected to predict the corresponding damage states of structures for socioeconomic loss computation [16].

2.2.2. Local Mechanism (Out-of-Plane Response)

Stone masonry buildings without rc slabs and ring beams having orthogonal walls not properly connected, or due to loss of in-plane walls integrity during seismic excitations, respond in local out-of-plane collapse of portion or complete walls [18]. The mechanical model for this mechanism is formulated elsewhere [21] for brick masonry walls (Figure 4).

In Figure 4, represents the peak acceleration of input excitation; represents the response acceleration of wall; represents total mass of wall; represents acceleration due to gravity; represents the wall thickness; , , and represent the limit state displacement capacities of wall.

The resistance of local out-of-plane mechanism to earthquake excitation is governed by the wall geometry, boundary condition, self-weight, and precompression level of rocking portion of wall while being less affected by the masonry material properties [22]. The out-of-plane stability of precracked wall can be relatively well assessed by using the secant stiffness, at the beginning of the third degrading branch of the mechanical model, and 5% damped response spectrum [22]. Thus, the secant period for out-of-plane mechanisms can be formulated as follows:where = represents the ratio of the limit state displacement capacity; = represents the collapse multiplier; represents the force at incipient rocking; represents the wall thickness; represents the total weight of the rocking wall. The ultimate displacement capacity of the out-of-plane wall is mainly governed by the thickness of the wall [21] which can be used with 5% viscous damping for seismic assessment and collapse prediction [22]. Nevertheless, the response of walls for out-of-plane loading in actual buildings is more complex than previously presented (Figure 4). Based on actual observations during earthquakes, various possible failure mechanisms are reported (Figure 5) [23] which are typical for most of the European historical towns. Collapse multiplier has been derived for different mechanisms using static analytical analysis to compute the lateral strength of out-of-plane loaded wall [23] which is further developed following detailed experimental investigation [24] on scaled models, typical for Italian dry stone masonry. The typologies of failure mechanisms for out-of-plane could be even more than the cases presented herein.

2.3. Development of Analytical Fragility Functions

Fragility functions, also called fragility curves, vulnerability curves, and/or building damage functions, describe the number of buildings reaching or exceeding a given damage state given the seismic intensity, represented as peak ground shaking parameters or spectral quantities [25]. Considering different possible damage states of buildings, there could be a number of fragility functions for a given typology which can be used to estimate the number of buildings in different damage level given the seismic intensity. Generally, for a given limit state, fragility function is derived considering a standard normal cumulative distribution function of the logarithmic difference of the seismic intensity and threshold capacity of limit states with certain level of standard deviation:where [..] represents the probability of reaching or exceeding a given limit state; represents the standard normal cumulative distribution function; SD represents the seismic intensity/demand; represents the limit state capacity of the system; represents the natural logarithmic standard deviation which defines the level of uncertainties in the fragility functions. The limit state capacity , usually median value, is obtained experimentally or numerically using sophisticated numerical tools. The standard deviation β is obtained from the square root square sum or similar combination of individual uncertainties which do not have clear rationale and justification behind. Similar other procedures exist which make use of constraint criterion to derive analytical fragility functions.

For loss estimation on regional scale, the uncertainties and variabilities in structural characteristics, geometrical and material uncertainties, can be obtained through on site survey of building stock and laboratory investigation of structural materials. The survey can better provide information on the likelihood of different geometrical features of regional building stock, for example, beam/column depth, width, length, reinforcement details, and number of structure’s storeys [16], which affect the seismic response of structural systems. No any hypothesis can give clear and conceptual information on how these local structural uncertainties can be truly represented in the fragility functional form, that is, , using (5). Thus, a transparent and conceptual methodology is presented to derive analytical fragility functions for regional building stock, taking into account different sources of uncertainties explicitly, without making any constraint.

2.3.1. Description of the Proposed Methodology

The DBELA methodology makes use of displacement-based analytical fragility functions where the seismic intensity is defined as a vector-based inelastic displacement demand on the structural system due to its direct correlation with the structural expected performance level for a given seismic demand [8, 20]. The seismic intensity is obtained in an iterative fashion following the damage states of structural systems for a given displacement response spectrum. The random linear 5% displacement response spectra are used to assess the structural performance and derive fragility functions which have the advantage to avoid the problem of being specific to a given regional spectrum and which can be used later for any risk assessment study using scenario earthquakes, uniform hazard spectra, and/or code spectra. To better understand different steps involved in the methodology flow charts are provided for fragility function’s derivation for both global and local mechanisms (Figures 6 and 7) while each of the major steps involved is described as follows.

(1) Global Mechanism

Step 1. The first step of the method is the generation of random population of buildings which represent a given class of building typologies within a given urban/rural exposure. Controlled Monte Carlo simulation is used to generate thousands of buildings, each with different geometrical and mechanical properties being defined using a complete probabilistic distribution with prescribed mean and coefficient of variation.

Step 2. The second step of the method is to define random seismic demand on the generated buildings which is performed through the use of random linear 5% damped displacement response spectrum. Special variability of ground motions is not considered in the fragility derivation which can be considered later in the application of fragility functions for developing damage scenarios for regional risk assessment and loss estimation.

Step 3. For each of the spectra, (a) for each limit state the secant vibration period, displacement capacity, and viscous damping of the buildings from the random populations are computed using calibrated structure-specific mathematical models, (b) the displacement demand on each of the buildings is obtained from the overdamped displacement response spectrum at the limit state vibration period of that building which is then compared with the displacement capacity of the building to predict its performance, and (c) the number of buildings having capacity less than the demand is summed and divided by the total number of the generated buildings to obtain the probability of exceedance of a given limit state. Similar hypotheses are used by other support methodologies [810] to compute the limit state probability of exceedance.

Step 4. For each of the spectra, (a) the median yield vibration period and median yield displacement capacity are obtained from the generated building stock, (b) the spectral displacement demand at the median yield vibration period is obtained from the 5% damped elastic displacement spectrum, (c) this is compared with the median yield displacement capacity, if in case the demand is less than the yield displacement (i.e., this defines the median spectral displacement demand on the structures), and (d) for spectral displacement demand greater than the median yield displacement, the performance point is obtained in an iterative fashion which defines the median spectral displacement demand (seismic intensity).

Step 5. The probability of exceedance for each limit state is plotted versus the median displacement demand (seismic intensity) for each of the random spectra. Available cumulative distribution functions are fit to the data and the unknowns of the functions are obtained to completely describe the fragility functions for future applications.

(2) Local Mechanism

Step 1. The first step of the method is the generation of random population of buildings which represent a given class of building typologies within a given urban/rural exposure. Controlled Monte Carlo simulation is used to generate thousands of buildings, each with different geometrical and mechanical properties being defined using a complete probabilistic distribution with prescribed mean and coefficient of variation.

Step 2. The second step of the method is to define random seismic demand on the generated buildings which is performed through the use of random linear 5% damped absolute displacement response spectra [20].

Step 3. For each of the absolute spectra, (a) vibration period, displacement capacity, and viscous damping of the buildings from the random populations can be computed for a given local mechanism at the collapse limit state using simplified empirical models, (b) the displacement demand on each of the buildings for local mechanism can be obtained from the 5% damped absolute displacement response spectrum at the global secant vibration period of building which is then amplified with the out-of-plane amplification spectrum/function at the vibration period of the considered local mechanism in order to compare the demand with the displacement capacity to predict its performance, and (c) the number of buildings having capacity less than the demand is summed and divided by the total number of generated buildings in order to obtain the probability of exceedance of local collapse mechanism.

Step 4. For each of the absolute spectra, (a) the yield vibration period for global response and median secant vibration period for local response are computed and (b) the spectral displacement demand at the median yield vibration period of the given class is obtained from the 5% damped absolute displacement spectrum which is amplified with the floor amplification function at the local mechanism’s secant period in order to obtain the median spectral displacement demand (seismic intensity).

Step 5. The probability of exceedance for each limit state is plotted versus the median displacement demand (seismic intensity) for each of the random spectra. Available cumulative distribution functions are fit to the data and the unknowns of the functions are obtained to completely describe the fragility functions for future applications.

3. Experimental Test on Case Study Rubble Stone Masonry

3.1. Model Description and Test Setup

The DBELA assessment approach needs the basic material properties of masonry and geometric detailing to develop structure-specific mechanical models (in-plane/out-of-plane) for the assessment of global and local capacity of masonry structural systems. The essential material properties, compressive strength, elastic moduli, tensile strength, and so forth, of considered rubble stone masonry cannot be obtained reliably at the section level due to the difficulties in performing tests on masonry prism, reproducing the true replica of the filed condition and huge scatter in the observed behavior, which is due to the fact that when stone-to-stone contact is found during the compression tests, very huge value of compressive strength is achieved which gets minimal value when stone-to-void possibility is found. Huge uncertainties in the material properties at the section level make the global response less reliable to be achieved. Also, the material properties at section level cannot directly provide guidance on the performance levels, damage scale, and predominant seismic response mechanism of a building. Thus, a full reduced scaled structural model is tested to obtain the global capacity parameters and develop damage scale in order to develop the DBELA tool for global seismic risk assessment of stone masonry buildings with rc slab in Pakistan.

One-third scaled structural model of single storey and single room is designed for corresponding prototype structure using the similitude principles (Figure 8). A comprehensive detail of the testing program and test results are published elsewhere [26]; the present paper will briefly describe the overview of the test and the main findings important within the scope of the present research work. The model is constructed of half-dressed stone masonry work in cement mortar with rc slab and ring beam which represents most of the existing urban building stock, public and private, in the region. The structural model is fixed to 130 mm thick concrete pad at the bottom which is mounted on the shake table, a single degree of freedom (SDOF) shake table at the Earthquake Engineering Center of University of Engineering and Technology Peshawar Pakistan, and firmly secured using bolts. The structural model is tested in the weaker direction using natural accelerogram, Kobe 1995, and incremental excitations with linear scaling until the complete collapse of structural model. The model is instrumented with accelerometers and displacement transducers at the base of the model and top of the roof: two accelerometers at the base and one on the roof, two displacement transducers at the base, and two on the roof right above the walls.

3.2. Observed Response of the Structural Model

The structural model is subjected to 5% scaled Kobe record with linear amplitude scaling until the total collapse of structural model. The model is inspected for the observed damage after every run. The floor acceleration and displacements, on each in-plane wall and at the base of the model, are recorded and processed. The floor acceleration is normalized by the seismic mass of the structural model in order to obtain the base shear for the corresponding equivalent SDOF system. Also, the average floor displacement demand is obtained which is corrected with the base displacement demand in order to obtain the relative displacement demand on the system, which is normalized by the model height to compute the corresponding drift demand on the system. The observed physical damage states of tested structural model are shown (Figure 9) at different run of input excitations and the capacity parameters, for the corresponding prototype buildings, along with the description of observed damage at different performance levels which are provided in Table 1.

4. Derivation of Mechanical Models for Stone Masonry Buildings

4.1. Characteristics of the Case Study Buildings

Stone masonry buildings are practiced abundantly in the northern areas, urban and rural, of Pakistan due to large local availability of stone material and low cost of labor in construction. Stone masonry with earthen floor, wooden/steel truss with GI sheet, and/or rc floor are the common residential building construction practice for single storey in rural areas and up to three storeys with rc floors in urban areas. Stone masonry in cement mortar contributes 50% in overall to the total building stock in the recent earthquake affected areas (http://www.erra.pk/). The most prevailing building dimensions range from 8 m × 5 m to 15 m × 5 m with typical wall density, with the ratio of the cross-sectional area of the in-plane walls to the total floor area, ranging from 10% to 15% [27]. The buildings have 300 to 500 mm thick load bearing walls with rubble stone masonry, having 130 to 150 mm thick rammed earthen roof, wooden/steel truss GI sheet, or rc slab and interstorey height of 2.0 m to 3.0 m. The buildings are provided with a ring beam right above the walls, approximately 150 to 230 mm deep, and width equal to the wall thickness. The building rests on shallow strip type footing, with stepped stone work overlain compacted earth surface. The load bearing walls are perforated by doors, with typical dimensions of 1.00 m × 2.13 m, and windows, with typical dimensions of 1.00 m × 1.22 m to 1.83 m × 1.22 m, thus providing 70% effective cross-sectional area of walls for carrying both gravity and lateral loads. The primary seismic resistance mechanism for stone masonry buildings with rc slab, as observed in the dynamic test, is in-plane global mechanism with diagonal shear cracks followed by bed-joint sliding in the short piers. However, the ultimate mechanism, collapse of the building, is governed by the combined in-plane and out-of-plane failure of masonry walls (Figure 9). Stone masonry buildings, which do not respect the minimum requirements to ensure in-plane integrity of walls, respond in local out-of-plane failure of portion of walls or complete walls and can be followed by the complete collapse of building [18].

4.2. Prototype 2D Cases Study Structural Models
4.2.1. Design of Prototype Buildings

Due to unavailability of basic material properties for stone masonry, and the reasons mentioned earlier, it is not straightforward to develop numerical tools for future applications for global analysis. Thus, the present study considered the simplified equivalent frame approach, SD-SAM [11, 12] earlier proposed and employed in the regional seismic loss estimation study, for nonlinear static and dynamic time history analysis (global analysis) of masonry buildings, and developed herein for case study buildings of Pakistan. The method uses the idea of modeling masonry spandrels and piers as equivalent beam-column elements with bending and shear deformation with infinitely stiff joint element offsets at the ends of the pier and spandrel elements (Figure 10), as proposed earlier [28]. The effective deformable length of pier is approximately obtained taking the intersection of the element with a line making a 30-degree angle with the corners of the openings. Each of the elements is provided with lumped nonlinear hinge having an appropriate force-displacement response to simulate the lateral response of masonry wall, depending on the elements ultimate mechanism. The properties of hinge, force-displacement response, are obtained from the available empirical rules and/or experimental tests on masonry material and walls [10, 11]. The present study performed the calibration of the SD-SAM approach for stone masonry with rc slab (global mechanism) as follows.

In the first step, a generic prototype of equivalent frame is generated for the tested structural model and designed, respecting the geometric detailing and loading conditions, with the available empirical shear strength models for masonry walls, previously developed [17, 29] and employed for brick masonry buildings [11, 12]. The beam-column element used for masonry idealization is completely defined by masonry Young modulus, shear modulus, wall sectional area, and the wall moment of inertia. The hinge is assigned with nonlinear force-displacement constitutive law to simulate the lateral capacity of the frame elements, that is, pier and spandrel. The total lateral displacement capacity of the element at any instant is the summation of the bending, which is contributed by the elastic flexure bending of beam-column element, and shear, which is contributed by the lateral translation of hinge, deformation developed in the element [11, 12]. The possible failure mechanisms considered for masonry piers are flexure mechanism and shear failure using the following strength models:where represents the ultimate strength for flexure failure mode; and represent the length and thickness of the wall; = represents the mean vertical stress due to axial load ; represents the total height of the pier; is 1.0 for a cantilever pier and 0.5 for a pier fixed-fixed boundary conditions; represents the compressive strength of masonry; represents the ultimate strength for diagonal shear failure mode; represents the diagonal tensile strength of masonry; for , for , and for ; represents the ultimate strength for sliding shear failure; , assumed as 0.4, and , taken as , represent the coefficient of friction and cohesion of masonry as global strength parameters. The strength computed is reduced by 10% in order to be conservative and respect the energy balance criterion [17, 29] which is observed to predict well the equivalent strength of masonry buildings [10].

In the second step, the tensile strength, , of masonry is selected, fixing the compressive strength to 4 Mpa, to achieve the equivalent yield strength, 4.50 m/sec2, of the corresponding prototype systems. In the next step, the Young modulus, , and shear modulus, , are selected to achieve the corresponding drift limits at yield and the 1st modal period of structural system. The assumption made herein for compressive strength and shear modulus is typical for such type of masonry system [30]. The equivalent frame is designed, using SD-SAM, employed in OpenSees [31] with masonry walls represented by equivalent elastic beam-column element which is provided with zero-length element, hinge (Figure 10), with hysteretic material to define the nonlinearity, force-displacement constitutive law, of masonry walls with bilinear rule. The structural model is analyzed to obtain the capacity curve with displacement controlled pushover analysis which is compared with the observed capacity of the tested prototype building and is found reasonable in terms of yield stiffness and equivalent bilinearised strength, respecting the energy balance criterion [17, 29] (Figure 11). The ultimate displacement capacity, and the corresponding drift limit of building, of the nonlinear hinge is obtained, when the lateral strength obtained from experimental results (see Table 1) drops by 20%.

4.3. Dynamic Analysis of Case Study Structural Models

The present study considered 125-case study 2D two-storey structural models, selected with different floor area, wall density, and material properties, in order to take into account the regional geometric and material uncertainties explicitly in the capacity evaluation, designed with the material sectional properties obtained in the previous section. The frame elements are assigned with bilinear Takeda type rule [32] (Figure 11) with Emori and Schnobrich [33] type of unloading stiffness and beta equal to 0.6 as proposed earlier [11, 12] for masonry buildings. The case study structural models are analyzed dynamically using nonlinear time history analysis (NLTHA) with 10 natural accelerograms extracted from the PEER NGA database for soft soil condition with the mean spectrum compatible to EC8 TypeI-C-soil spectrum [34] (Figure 12).

The accelerograms are linearly scaled in order to observe the postyield response of the models which is used then to derive static SDOF system, mechanical models, for the considered building typology in order to retrieve the dynamic characteristics of buildings in preyield and postyield limit states and develop empirical models for limit state’s vibration period and displacement capacity. It is worth mentioning that the numerical tool used herein is simplified but respecting the fundamental structure’s parameters, that is, stiffness, strength, ductility, and energy dissipation, sufficient to demonstrate the global structural performance for seismic loading. Also, the findings of numerical tool are used for the risk assessment of buildings on regional scale where, apart from large case study structural models, the regional structural geometric and mechanical properties are approximately defined which makes the present tool as a reasonable choice for structural analysis.

The scope of the dynamic analysis is to compute the equivalent base shear and equivalent displacement demand at different limit states of masonry walls on the ground floor, since ground storey mechanism is governed for the considered typology, using the proposed SDOF derivation [11, 12], and computes the secant vibration period of considered building stock:where represents the equivalent displacement of structural system; represents the equivalent base shear for the corresponding SDOF system; and represent the floor mass and floor displacement demand, respectively; VB represents the base shear demand; represents the equivalent mass of the structural system. The vibration periods obtained for the considered case study buildings are used to develop the period model (10) for future applications in regional risk assessment using the DBELA methodology. where represents the total height of the building; and are the coefficients obtained through regression analysis of the data; represents the regional variability in the period computation due to geometric and material uncertainties and record-to-record variability. The period coefficient with set to 0.75 obtained for 25 structural models with record-to-record variability and mean material properties is shown in Figure 13. Also, the period coefficients obtained for all the 125 structural models are obtained.

It is observed that for a given structural system the period reduces with increasing wall density. Similar trend of vibration period is also observed for increasing floor area. The period coefficient obtained for all the cases study structural models, 1250 cases, takes into account the geometric (floor area and wall density) and material uncertainties besides the variability introduced by earthquake excitations. Existing and conventional code-based formulae considered in (10) equal to zero which gives the mean estimate of vibration period. However, the height dependent only formulae will neglect the floor area and wall density effect in the buildings period estimation for regional risk assessment.

Additionally, the structural models are analyzed to develop the limit states displacement capacity model (2) for different performance levels of the structures. In the first step, static analyses are performed for all 125 structural models to obtain the crack and yield limit state drift values using the shear strength, obtained using calibrated strength models ((6)~(8)), and stiffness at cracking and yielding, respectively, with 50% cracked sectional properties and 70% effective shear area, of masonry walls, respectively. The crack force is conservatively considered being 60% and the yield force being 90% of the maximum computed force [17, 29]. A mean value of 0.24% is obtained for crack limit state with logarithmic standard deviation of 0.16, min 0.16%, and max 0.33%. Similarly, for yield limit states, a mean value of 0.35% is obtained with logarithmic standard deviation of 0.15, min 0.24%, and max 0.50%. These values are less than the experimental results due to the fact that the tested model has wall density of 12.84%. The computed drift values take into account the geometric and material variability for rather broader range of wall density and floor area. Using the limit states ductility obtained in the experiment, which is reduced by 20% in order to be conservative, a mean value of 0.768% and 1.016% is obtained for major and collapse limit states, respectively, which nevertheless need to be updated once more experimental data on rubble masonry are made available. Also the structural models are analyzed to obtain the displacement coefficients, and , at the preyield and postyield limit states, which are used herein to take into account the record-to-record variability in the displacement capacity evaluation of the SDOF systems and make the model (2) treat different uncertainties, that is, geometric, material, and record-to-record, explicitly in the capacity evaluation. Each of the structural models is investigated to obtain the floor displacements and the corresponding base shear at the yielding of the ground floor walls for each NLTHA of a given structural model which is used to obtain the effective height of the SDOF system, following the recommendation of Priestley et al. [20]:where represents the effective height of the SDOF system, which corresponds to the center of seismic force; represents the floor height. The displacement coefficient at yield limit state () is obtained then by dividing over the total height of the structural model, . A mean value of 0.7269 is obtained with logarithmic standard deviation of 0.0115, and minimum value of 0.6880 and a maximum value of 0.7500 are obtained. Since the considered building typology has soft-storey ultimate mechanism with demand mostly concentrated at the ground storey level, thus the above hypotheses cannot be directly used to obtain displacement coefficient at postyield limit states (); hence the following procedure is used. NLTHA of all the structural models is performed, through linear scaling of accelerograms, to exceed the near collapse limit state of walls at the ground floor. The data is analyzed to obtain the equivalent displacement demand when the target ductility of 3 is exceeded at the ground floor level, which corresponds to the collapse limit state of structural model. In the next step (2) is used, inverting and rearranging, to obtain :where represents the equivalent displacement at the near collapse limit state of ground floor walls; represents the equivalent displacement at the yielding limit state of ground floor walls; represents the plastic drift demand on the ground floor walls, obtained by dividing the ground floor displacement demand over the ground floor height; represents the ground floor height. A mean value of 0.8068 is obtained with logarithmic stand deviation of 0.1415, and minimum value of 0.439 and a maximum value of 1.2335 are obtained. Figure 14 shows the limit state displacement profiles of a two-storey structural model, 40 m2 floor area and 10% wall density, for 10 accelerograms and the displacement coefficients, and , for all the 125-case study structural models.

5. Derivation of Analytical Fragility Functions

The proposed methodology is used herein to derive analytical fragility functions for urban rubble stone masonry buildings of Pakistan. Application of the methodology is also carried out for brick masonry urban building stock of NE Pakistan [35]. However application on the derivation of fragility functions for rural stone masonry building stock is shown also herein which exhibits mainly out-of-plane local failure of walls. The following section will briefly describe the global and local fragility function’s derivation.

5.1. Displacement-Based Analytical Fragility Functions for Stone Masonry Buildings of Pakistan
5.1.1. Global Mechanism

Controlled Monte Carlo simulation is used to generate random building stock with different geometric and mechanical properties considering lognormal probability density function (pdf) for all the parameters involved in the capacity evaluation (Table 2). The lognormal pdf is considered for simplicity reasons and to be conservative in structural capacity estimation. The damage scale is selected as observed in the experimental test, Table 1. The displacement capacity and secant vibration period are obtained at each limit state, using the calibrated empirical models, that is, (1)~(2). Due to unavailability of material viscous damping, the damping model proposed for other masonry types [11] is considered herein.

5.1.2. Local Mechanism

Stone masonry buildings with rc slab did not show prominent out-of-plane failure in experimental test and during 2005 Kashmir earthquake but were damaged mainly due to diagonal cracks in wall, cracking around the corners of opening, horizontal cracks at the roof level, and complete collapse of buildings [36]. However stone masonry buildings with wooden/steel truss and GI sheet suffered major local out-of-plane failure of portion of walls or complete walls. Collapse of gables, corner wedge type of wall failure, and complete overturning of exterior walls were the common out-of-plane failures [36]. Delamination of exterior wythe of walls was observed in few cases but this never caused serious problem to the overall stability of buildings [36]. Two-way vertical bending [22] of out-of-plane walls is not observed; thus the mechanical model (4) derived for the free standing out-of-plane rocking walls is considered only. The limit state displacement capacity of out-of-plane mechanical model is adopted from Doherty et al. [21]; however the ultimate displacement capacity of wall is considered with reduction factor () in order to be conservative in capacity evaluation [8, 24, 37]. Following experimental investigation on the out-of-plane testing of several dry stone masonry walls at the University of Pavia [37], a reduction factor () is proposed for the collapse multiplier obtained using rigid block behavior of out-of-plane responding walls. Different parameters involved in the in-plane/out-of-plane mechanical models are defined to generate regional building stock, both urban and rural building stock. Table 1 shows the values assigned to each parameter and the assumed distribution functions to generate random numbers using controlled Monte Carlo simulation.

Once the regional building stock is generated and the limit state capacities are evaluated in a probabilistic fashion, random linear 5% damped displacement response spectra are generated. The global in-plane assessment is performed through overdamping the linear spectrum, using overdamping factor [38] and the system viscous damping, for each limit state. where represents the overdamping factor; represents the masonry viscous damping [11]. The local out-of-plane assessment is performed using absolute displacement spectrum along with the out-of-plane amplification spectrum [20]:where represents the period of out-of-plane wall; represents the in-plane secant period of building; represents the viscous damping of out-of-plane wall, recommended as 5% [22]. The relative and absolute displacement spectrum and out-of-plane amplification spectrum used for in-plane and out-of-plane assessment are depicted in Figure 15.

The absolute spectrum is anchored at the peak ground displacement (pgd) and linearly increased to maximum displacement demand at corner period following the recommendation of Priestley et al. [20]. Different recommendations can be used to compute the pgd, corner period, and corner displacement for worldwide regions [34, 3944]. However, for the out-of-plane fragility function’s derivation the pgd and the slope of absolute displacement spectrum are randomly selected using Monte Carlo simulation while increasing the spectrum linearly, against the period, without any bound for convenience and to be conservative. Nevertheless, an appropriate displacement spectrum, both relative and absolute, has to be used for developing damage scenarios for intended loss estimation study. The procedure outlined in Figures 6 and 7 is used to derive the fragility functions for two-storey stone masonry building stock of Pakistan for both in-plane global mechanism and local out-of-plane mechanism as shown in Figure 16.

5.2. Structure-Independent Fragility Functions

It is mentioned that fragility functions are needed to be developed in terms of structure-dependent vector-based intensity measures (inelastic displacement demand in the present case which is well correlated with the structural expected performance) and in order to avoid the problem of being specific to given assumed spectral shape and region, for example, code spectra and UHS for a given site which are not conceptual [45, 46] when losses are required for scenario earthquakes. However, for decision-making it will be of utmost importance to present the earthquake impacts as a function of convenient and easy to perceive structure-independent parameters, yet respecting the fundamentals of strong ground motions (frequency content, ground motion uncertainties, site soil characteristics, and correct spectral shape) and structures (strength, stiffness, ductility, and energy dissipation). Also, a prompt assessment for quick response and emergency planning soon after the earthquake could be much facilitated using fragility functions which need fewer parameters to estimate the socioeconomic impacts of an earthquake event. The simplest and readily available parameters to compute earthquake losses at a given site for an earthquake event are the magnitude of earthquake, source-to-site distance, site soil condition, site building typology, building content, and occupancy.

Thus fragility functions are derived for the MDR of stone masonry buildings of Pakistan expected at a site for a given scenario earthquake. Different possibilities of source-to-site distance and magnitude are selected, considering soft soil condition (type D of NEHRP soil classification as recommended for Pakistan [47]). The empirical ground motion prediction equation [44] of PEER NGA project [48] is used to generate ground motions for a given scenario earthquake with due consideration of ground motion uncertainties (both inter- and intraevent), 10,000 simulations for each scenario earthquake. The selected GMPE has the advantage to predict also pgd for scenario earthquakes. Figure 17 reports the structure-independent fragility functions for the considered building typology of Pakistan. The earthquake loss model of Bal et al. [16] is used to compute MDR. For a given site, knowing the MDR, the unit replacement cost of a building, and the total number of buildings in a given municipality, the expected economic loss can be easily obtained [11]. Different other loss models are available to compute the expected human causalities in heavy damaged and collapsed buildings [49]. However, fragility functions are derived herein for the trapped people along with their expected level of injuries using site-specific casualty model, derived from the recent 2005 Kashmir Earthquake [50].

As expected from the structural mechanics standpoint, the local out-of-plane mechanism results in relatively higher seismic risk and losses as compared to global in-plane mechanism of structural systems. Depending on the site building characteristics, if buildings are provided with rc slab and ring beams (like the present situation in the earthquake affected areas, http://www.erra.pk/, and some urban exposure) only the global mechanism’s charts have to be used to compute the regional seismic risk and losses for a given scenario earthquake, described only by the source-to-site distance and moment magnitude. If the regional buildings do not meet the minimum criterion to ensure the in-plane integrity of structural system and global seismic response mechanism [18], only the out-of-plane mechanism’s charts have to be used for assessment. A weighting factor has to be applied to MDR in cases where both in-plane and out-of-plane mechanisms are expected in regions; that is, for a given scenario MDR can be obtained for both mechanisms which will be multiplied by the corresponding percentage of each buildings class (global mechanism/local mechanism).

The expected causality rate is presented in terms of the occupancy in a housing unit and the total number of building stocks in the region. Thus the number of housing units in a region and the likelihood function of occupancy per housing unit, 7 to 12 for the considered typology [27, 51], have to be multiplied with the number obtained from the casualty rate charts (Figures 18 and 19).

5.3. Comparison with the Observed Response during 2005 Kashmir Earthquake

The structure-independent fragility functions are used to assess the seismic risk of stone masonry buildings for a scenario earthquake that recently occurred in the country, that is, 2005 Kashmir Earthquake. The Kashmir earthquake has a moment magnitude of 7.6 that struck most of the northern areas of Pakistan on 8 October 2005 at 08:50 AM local time (https://earthquake.usgs.gov/earthquakes/eventpage/usp000e12e#executive). The epicenter of the main earthquake shock was located at Latitude 34.49N and Longitude 73.63E, approximately 19 km north east of Muzaffarabad city of Azad Jammu Kashmir. The earthquake event is the direct results of collision of the Indian plate, moving 40–50 mm/yr northward, with the Eurasian plate. The collision is crust to crust and thus produces large magnitude shallow earthquakes [52]. More than 780,000 buildings were either destroyed or damaged beyond repair during Kashmir earthquake and more rendered unusable and needed demolishment and replacement (http://www.eeri.org/lfe/pdf/kashmir_eeri_2nd_report.pdf). The destruction was observed over an area of 30,000 sq-km causing severe loss of life and property [53].

Almost all the buildings, mainly stone and block masonry laid in cement sand mortar with rc slabs or GI sheet roofing, collapsed in the areas close to the epicenter [2]. In regions approximately 25 km away from the epicenter nearly 25 percent of the buildings collapsed and 50 percent of the buildings were severely damaged which gives a MDR of about 0.79 using the loss model of Bal et al. [16]. To test the proposed methodology and structure-independent scenario-based fragility functions, seismic risk is computed for considered buildings considering both local and global vulnerabilities. The structure-independent fragility functions are derived for individual mechanisms and also for the combined mechanism (Figure 20). The combined fragility function is derived using the weighting factor 0.20 for in-plane mechanism and 0.80 for out-of-plane mechanism as the building stock proportion in the earthquake affected regions was found to be 20 percent urban buildings (represented by in-plane mechanical model) and 80 percent rural buildings (represented by out-of-plane mechanical model) [3]. The site soil is considered as type D of NEHRP soil classification with shear wave velocity of 250 m/sec2.

Using only the in-plane fragility function gives an estimate of the MDR equal to 0.26 which gives an error of about +67.09 percent thus underestimating the risk, and the out-of-plane fragility function gives estimate of MDR equal to 1.00 which gives an error of about −26.58 percent thus overestimating the seismic risk, and the combined fragility function gives an estimate of MDR equal to 0.85 which gives an error of about −7.60 thus overestimating the site risk. Nevertheless, the observed MDR is reasonably predicted in case of combined fragility function. Also, the proposed fragility functions are conservative, that is, giving slightly overestimated seismic risk and losses.

5.4. Uncertainties in MDR

The mean prediction for seismic risk can reasonably predict the observed response during earthquake, as observed in the previous section. Nevertheless, considerable uncertainties can result in loss estimation for a given region considering a single scenario earthquake [54] and/or fragility functions [55]. The present methodology is capable of providing analytical estimate of the uncertainties expected in the regional seismic risk assessment. Thus the uncertainties are quantified in MDR to provide better guidance in decision-making. Each of the scenarios is simulated considering 10,000 ground motion possibilities with the uncertainties defined by the ground motion prediction equation (total uncertainties are used) [44]. The coefficient of variation (COV) and standard deviation are obtained for each scenario and plotted to observe the trends, and quantify, in COV and standard deviation against the MDR (Figure 21).

It is observed that COV tends to decrease with increase in MDR, a similar trend also observed elsewhere [55], and there is little dependence on the magnitude of earthquakes. On the other hand the standard deviation is found to behave differently where initially the standard deviation increases with increasing MDR, as observed elsewhere [55], up to certain level but then tends to saturate when the MDR gets very higher. The values of COV and standard deviation at a given MDR are different than the reported one [55]. It is worth mentioning that other reported findings [55] are for different regions and structural systems which gives an indication that COV and standard deviation in loss estimation studies are structural and regional dependent.

6. Conclusions and Future Developments

Seismic risk assessment of stone masonry buildings of Pakistan is performed using an innovative and state-of-the-art nonlinear static methodology for seismic risk assessment and loss estimation of buildings on regional scale. The methodology takes different possible sources of uncertainties in an explicit and transparent manner to compute the structural capacity and assess the seismic performance of structural system. Analytical displacement-based fragility functions are derived for the considered buildings, for both local and global vulnerabilities, which are used to estimate the seismic risk and losses for scenario earthquakes and develop structure-independent fragility functions. Simplified charts are provided to compute the socioeconomic impacts of earthquakes given the magnitude and source-to-site distance of the event. Comparison of the methodology and derived fragility functions is performed with the recent earthquake observation, in predicting seismic risk of the considered typology, which is reasonably predicted. The methodology can provide also estimate of uncertainties in regional seismic risk to provide help in decision-making. However, from the derived structure-independent fragility functions it is possible to develop seismic risk prediction equations for regional loss estimation. Comparison with large earthquake observation databases can provide estimate of inter- and intraevent uncertainties in seismic risk prediction in similar fashion as ground motion prediction equations do. It will make the regional loss estimation studies to assess the regional risk directly using the derived equations, once developed for a region. Similar studies will be performed on other regional building stocks of Pakistan once more experimental and reliable data are made available for the region. Additional data on stone masonry material and buildings can further improve the findings herein in order to increase confidence in future applications. The findings herein are the first of its kind for the considered region.

Disclosure

The research work presented herein is the further extension of the methodology, adopted for the fragility functions derivation under the EU funded project SYNER-G (Systematic Seismic Vulnerability and Risk Analysis for Buildings, Lifelines Networks and Infrastructures Safety Gain), to derive functions for the direct seismic risk assessment of buildings given the scenario earthquake.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.