Abstract

A stochastic-event probabilistic seismic hazard model, which can be used further for estimates of seismic loss and seismic risk analysis, has been developed for the territory of Yemen. An updated composite earthquake catalogue has been compiled using the databases from two basic sources and several research publications. The spatial distribution of earthquakes from the catalogue was used to define and characterize the regional earthquake source zones for Yemen. To capture all possible scenarios in the seismic hazard model, a stochastic event set has been created consisting of 15,986 events generated from 1,583 fault segments in the delineated seismic source zones. Distribution of horizontal peak ground acceleration (PGA) was calculated for all stochastic events considering epistemic uncertainty in ground-motion modeling using three suitable ground motion-prediction relationships, which were applied with equal weight. The probabilistic seismic hazard maps were created showing PGA and MSK seismic intensity at 10% and 50% probability of exceedance in 50 years, considering local soil site conditions. The resulting PGA for 10% probability of exceedance in 50 years (return period 475 years) ranges from 0.2 g to 0.3 g in western Yemen and generally is less than 0.05 g across central and eastern Yemen. The largest contributors to Yemen’s seismic hazard are the events from the West Arabian Shield seismic zone.

1. Introduction

The tectonic movement and interaction of the Arabian and the African plates constituting rifts of the Red Sea and Gulf of Aden are the principal cause of earthquakes in Yemen, which are known from historical sources over the last millenniums [13]. Earthquakes that affect Yemen are mainly associated with rifts of the Red Sea and the Gulf of Aden [4]. However, small- to moderate-sized earthquakes, which occur inside the Arabian Plate within 200 to 300 km of the axis of the Red Sea [5, 6], may also have some impact on the territory (Figure 1).

In general, seismic history of Yemen indicates the occurrence of large earthquakes with 20- to 30-year recurrence periods [2]. We can mention the Dhamar earthquake of 13 December 1982 (), which was felt over a large area and which killed and injured more than 15,000 people and destroyed about 1,500 settlements, as the deadliest earthquake of Yemen [8]. Two significant earthquakes occurred in 1941 and 1909. The main event () and two aftershocks ( and ) in 1941 killed about 1,200 people and damaged about 1,400 homes. In the 1909 earthquake event, about 300 people were killed and approximately 400 houses were destroyed.

Seismic hazard analysis for the territory of Western Arabia and Yemen was performed by several researchers [911]. The conventional Cornell approach [12] has been applied in almost all studies, and the results were not suitable for seismic loss analysis. The goal of this study is the development and testing of a stochastic-event probabilistic seismic hazard model for the territory of Yemen, which considers the influence of local site conditions and which can be used further for estimates of seismic loss and seismic risk analysis.

2. The Approach

The information related to the expected seismic effect and expressed in terms of earthquake ground motion parameters, such as seismic intensity and peak amplitudes of ground motion, is necessary for estimation of seismic losses and design of buildings and structures in earthquake-prone regions. Thus, the specification of engineering ground-motion parameters is the goal of seismic hazard analysis.

Probabilistic seismic hazard analysis (PSHA) involves the quantitative estimation of ground shaking hazard at a particular site taking into account characteristics of potentially dangerous earthquakes around the site. The principles of PSHA are well established within the scientific community. In this study, the source zone approach together with stochastic event generation has been used for developing a probabilistic seismic hazard model for Yemen. The main features of the approach are as follows: seismic source zones, which are considered as a potential source of earthquakes in the future, are divided into a series of linear sources (faults) and further into a system of elementary segments. Possible earthquakes may occur along the segments, and earthquake occurrence is assumed to be a stationary random process. Every elementary linear segment is characterized by the following parameters: (a) magnitudes of individual earthquake event and (b) annual rate of earthquake recurrence. Parameters of earthquake ground motion are calculated for every stochastic event (earthquake generated by elementary segment) along the territory of interest. The scheme allows using different parameters of seismicity and earthquake characteristics for different zones, as well as various ground-motion prediction equations.

3. Earthquake Catalogue

A primary component of the PSHA model is the earthquake catalogue for the region. Yemen has a long and relatively well-documented history of important past earthquake events including historical chronicles [2]. As historical data is not based on instrumental observations, the earthquake magnitude, intensity, and epicenter in such database have been estimated based on the description in historical documents. The development in instrumental seismology in Yemen started after the occurrence of the December 13, 1982, Dhamar earthquake [3].

A composite and updated catalogue has been compiled for the present study using two basic sources. The first source is the database collected by Ambraseys et al. [7]. The data has been divided by the authors into two subsets, which cover two periods from 184 BC to 1899 AD and from 1899 to 1987. In the first subset, epicenters of earthquakes were determined from macroseismic data, and, therefore, in most cases they reflect centres of the most heavily affected area. The accuracy of earthquake location in such cases depends on the availability of observations. The second subset combines instrumental data with macroseismic information. The original locations of all events were reexamined and, in most cases, parameters of earthquakes were reassessed using teleseismic and macroseismic information. The second source of earthquake data is the World-ANSS (Advanced National Seismic System) worldwide catalogue hosted by the Northern California Earthquake Data Center. The ANSS catalogue has been created by merging the master earthquake catalogues from contributing ANSS institutions and by removing duplicate solutions for the same event. Besides these two main sources of the data, parameters of important historical earthquakes of Yemen were checked using other literature sources [2, 1315].

When compiling the composite catalogue in this study, the following has been considered. The catalogue of Ambraseys et al. [7] is a more authentic and reliable source of data than the other sources, because each event in the catalogue has been carefully reexamined. However, location and magnitude of many earthquakes depend entirely on macroseismic data and the event coverage is not uniform in time and area. Many events with no felt effects may have been neglected in the database. From 1963 onward, the ANSS catalogue has been compiled using the data from automatic recordings seismographs. Thus, the composite catalogue (Figure 2) has been created considering (1) all records from Ambraseys et al. [7] from 182 BC to 1962 and (2) all data from the ANSS catalogue from 1963 onward. Aftershocks were removed from the composite catalogue using the declustering method of Gardner and Knopoff [16]. The procedure defines a space and time window after each event, which is supposed to be a main shock [16, 17]. All subsequent events within this window are declared dependent events and omitted from the declustered catalogue—unless their magnitude exceeds the main-shock magnitude. Window thresholds and usually depend on the main-shock magnitude and they are assumed universal for the entire study region and study period.

Several magnitude scales (local magnitude , body wave magnitude , surface wave magnitude , moment magnitude , etc.) have been used by various researchers to characterize earthquakes in the region [18, 19]. In our study, we have considered moment magnitude , which is consistent with ground motion computation. Thus, if necessary, magnitudes in the composite catalogue were recalculated. The empirical equations suggested by Johnston et al. [20] which were developed using worldwide datasets, were used for and conversion; the empirical equation suggested by Ekström and Dziewonski [21] has been used for conversion (Figure 3).

The most important requirement for any earthquake catalogue is to judge completeness of the catalogue with respect to magnitude and time of event occurrence. Magnitude of completeness () is the lowest magnitude, at which 100% of the events are detected in a space and time. Below , some parts of events are missed in the catalogue due to the influence of several factors including the network limitation. A traditional method of frequency-magnitude distribution, which describes the relationship between the frequency of occurrence and magnitude of earthquakes, has been used to estimate , which describes the relationship between the frequency of occurrence and magnitude of earthquakes. The frequency-magnitude curves in Figure 4 show a higher frequency of earthquakes at certain magnitudes and after that it is decreased exponentially for higher magnitude and is picked up as 4.8 ± 0.1, at the point where the curve becomes exponential.

Although historical data goes as far back as 184 BC, the reporting is not homogeneous for the entire span of the catalogue to derive correct parameters. To compute catalogue completeness with respect to time and magnitude, events have been plotted for various broad regions at 0.2 and 0.5 magnitude interval ranges. It is observed that despite compiling the Yemen seismic regions catalogue from various sources, the earthquake catalogue is not complete uniformly for all magnitude ranges. Catalogue completeness has been computed following the method of Stepp [22]. For magnitude 6.0 and above, the catalogue is complete for the last 110 years, while for magnitudes between 3.5 and 4.0 the catalogue is observed to be complete for only the last 20 years. Table 1 shows the catalogue completeness times for Yemen seismic regions, which have been taken into account while computing recurrence parameters.

4. Seismic Sources

One of the main features of seismic hazard analysis involves identification of seismic sources. Seismic sources are geographical areas that have experienced seismic activity in the past and serve as potential sources of earthquakes in the future. Seismic sources are delineated based on seismotectonic features of the region and it is assumed that the past earthquake activity is a reliable predictor of the future activity.

A cluster of earthquakes occurred in the southwestern part of the Arabian plate near a complex triple junction of active spreading ridges along the Red Sea, Gulf of Aden, and Afar depression (Figure 1). In general, the seismic activity around Yemen is most pronounced along the spreading ridges. However, some small- to moderate-sized events occurred also within the Arabian plate at distances up to 200–300 km from the axis of the Red Sea [5].

A total of 12 seismic sources reflecting peculiarities of seismicity have been delineated in the studied area on the basis of seismotectonic maps, clustering of epicenters, and smoothening of seismicity (Figure 5). Also, the results of the Global Seismic Hazard Assessment Programme [23] and corresponding studies [24] were considered. Due to errors associated with locations of preinstrumental earthquake events and catalogue incompleteness, both the seismotectonic province approach and the epicenter smoothing Frankel approach [25] have been used to demarcate the area of source zones.

The epicenter-smoothing approach [25] has been considered as a lower bound estimator for seismic hazard and helps in decision making in moderate seismicity regions where source zone definition and estimation of maximum possible magnitudes can lead to a wide variety of estimates. In each broad seismic region, using the Gutenberg-Richter magnitude frequency distribution, -value (total number of earthquakes per year greater than ) and the slope (-value) at each 10 km grid cell have been computed. Then, using the Frankel approach, the cumulative annual earthquake occurrence rates have been computed in 0.1 magnitude increments in each 10 km grid with a correlation distance of 50 km [25, 26]. A thematic map of the smoothed -value at the grid level is presented in Figure 5 together with the defined seismic sources.

In each seismic region, the parameters of the Gutenberg-Richter magnitude frequency distribution where is the cumulative number of events greater than magnitude and -value () and -value () are the regression parameters, were estimated.

For each source, the constants -value and -value of the recurrence relationship have been obtained from regression analysis on the yearly rate of occurrence of historical events (Figure 6). While estimating yearly rate of occurrences, catalogue completeness by magnitude bins has been considered. The rate of occurrence of future events is calculated by binning events into half magnitude increments. The rate of occurrence of events in a given bin has been equal to the difference between the cumulative rate of occurrence of that bin and the cumulative rate of occurrence of the next bin ( larger). The rate of occurrence of max has been calculated in the same manner, taking the next bin as max plus .

Demarcated seismic sources have unique tectonic regimes and are different from adjacent sources in terms of stress accumulations and event recurrence characteristics. Based on clustering and homogeneity in historic seismic activity, earthquake source zones that have similar geophysical characteristics within a geographic area have been delineated. Ideally, each area source should have been considered as an independent seismic source to compute recurrence parameters. However, difficulties were encountered due to low seismicity for some sources. Some seismic areas especially within the Yemen Arabian Plate have sparse distribution of historical seismicity. The western region of Yemen has a high potential of seismic hazard as compared to the eastern part. In order to overcome the problem of low seismicity in the source area, events from certain sources have been merged for the purpose of estimating recurrence parameters. For example, parameters for seismic source Red Sea Al Darb have been computed by taking into account the individual source events as it had enough events for regression. However, five sources in the Arabian shield regions comprising sources 1 to 5 have been merged to derive the hazard parameters.

In estimating the seismic rate, we considered truncated exponential magnitude recurrence relationships (Figure 6) in which recurrence frequencies for upper magnitude are truncated [27]. This approach is important in order to (a) provide compatibility with the physical limitations of magnitudes that are specific to various seismogenic zones and (b) lead to realistic results of recurrence characteristics of local hazard at various locations, in cases when local hazard is analyzed (as it should be) on the basis of probabilistic convolutions between the magnitude recurrence laws and the random attenuation laws.

After recurrence relations were determined for the major sources with merged data the seismicity represented by relation was redistributed to each area source. In the process, it was assumed that the -value remains the same in all major seismic regions having association with 3–5 area-sources and only appropriate -value had to be assigned. Characteristics of the seismic zones are listed in Table 2.

5. Stochastic Events

Seismic area source zones have been used to model the earthquake occurrence in areas where observed seismicity exhibits a diffuse pattern and where specific earthquake-generating faults cannot be identified to link with historical seismicity. Thus, it is assumed that earthquakes can occur by rupturing any part of the source area zone. In performing seismic hazard analysis, normally many researchers including computer code such as EQRISK [28] are using point source for generating earthquake events by subdividing the area seismic zone into grids of equal sizes [29]. In reality occurrence of an earthquake event is never a point source, rather the source is always a fault rupture. This can be observed by elongated isoseismal lines and elongated pattern of aftershocks around an event. Instead of considering point sources, we have considered line sources (fault segments) in seismic hazard simulation. In our approach, the area zones have been subdivided into a series of line sources (fault segments) and possible seismic events occur along the line sources (Figure 7). Ideally, the length of the fault segments depends on magnitude-rupture length relationship and orientation should be along the identified natural faults. Since there is no comprehensive model available for the regions, in this study we assumed the same length for all segments. In reality, there is always a bias on distribution of natural faults within an area source. To avoid uneven seismic distributions, a systematic fault line along the preferred tectonic orientation has been assumed for stochastic events (Figure 7). Although the distribution and orientation of stochastic faults look unnatural, this assumption has two advantages: (a) occurrence of an earthquake event can be considered by rupturing a fault line, not a point source, and (b) seismic annual rate of occurrences can be distributed uniformly.

With an assumption that there is an equal probability for an event of particular magnitude to occur by rupturing any part of the delineated seismic area source zone, a number of stochastic faults have been compiled. Following the Poisson process, seismic activity of the given seismic source zone has been redistributed to the line source taking into consideration the length of the fault segment. These demarcated stochastic faults are treated as individual area line sources. The stochastic earthquake events set along these sources, with associated rate of occurrences, have been generated for earthquake simulation for ground motion footprints.

It is assumed that every segment can produce earthquakes within a considered range of magnitudes (from minimum considered magnitude min to maximum possible magnitude max). To capture possible scenarios, the set of stochastic events, which is associated with each line source, has been created at 0.2 magnitude interval starting from min = 4.5. A total of 15,983 stochastic events with their annual rate of occurrences, associated with 12 seismic sources and 1,583 fault segments, have been generated for hazard simulations. Table 2 shows the characteristics of seismicity assigned to seismic source zones and the corresponding number of segments and stochastic events.

6. Ground Motion Model

The level of ground shaking in seismic hazard analysis is usually expressed in terms of peak ground amplitudes (acceleration PGA, velocity PGV, or displacement PGD), spectral acceleration at discrete frequencies, and seismic intensities (on MMI or MSK intensity scale). The estimations have been obtained using an appropriate ground-motion prediction equation taking into account dependence between the given ground-motion parameter, earthquake magnitude , source-to-site distance , and local site conditions.

The unavailability of strong motion models for the Yemen region makes it necessary to select attenuation equations derived from statistically significant data sets from the literature. The collections of ground-motion prediction models, which were developed during the last decades, are described, for example, by Douglas [30, 31] and Abrahamson et al. [32]. However, the applicability of the generalized worldwide models should be carefully studied for the regions lacking the data to avoid additional uncertainty.

Ground-motion models in terms of PGA were developed for the region by several authors [33, 34]. The most recent model for the southern Dead Sea Transform region has been proposed by Al-Qaryouti [35] for the larger value from the two horizontal components as follows: where PGA is expressed in units of g and   is the closest distance to rupture. Thenhaus et al. [36] suggested a region-specific adjustment to the Campbell [37] relation to make it suitable for Saudi Arabia as follows: where PGA is expressed in units of g and is the closest distance to rupture. The model has been used by Al-Haddad et al. [24] for evaluation of seismic hazard and design criteria for Saudi Arabia. In developing seismic hazard and design criteria for Saudi Arabia Al-Haddad et al., 1994 had used these attenuation equations for all seismic sources in the vicinity of the Arabian Peninsula including Yemen. Youngs et al. [38] developed ground motion models for subduction interface zones and intraslab regions using worldwide data. Here, we used the following model that relates to subduction intraslab events and rock sites: where PGA (geometrical mean of two horizontal components) is expressed in units of g and is the closest distance to rupture, is the source depth, and SD is the magnitude-dependent standard deviation, .

Seismicity of Yemen is attributed mainly to earthquakes from the transform fault region and the stable continent region. Attenuation 1 is the most recent model for transform region, attenuation 2 is an established equation for stable Arabian continent, and (4) is the global model. In this study we selected these regional attenuation models, which represent ground motion for rock condition (2)–(4), to be applied in seismic hazard analysis for Yemen. In order to capture epistemic uncertainty in ground motion estimations, these three models were used with equal weight. The expected values of ground motion parameters from all stochastic events were estimated from averaged attenuation equations as shown in Figure 8.

7. Soil Modifications

It is well understood that near-surface geological conditions may strongly affect earthquake ground motion at a particular site amplifying the shaking amplitude and changing frequency of the motion. In certain cases (e.g., building code provisions), a few generalized site classes have been selected to describe the variety of local soil conditions and to characterize their effect on ground motion. A widely used site classification system has been based on the properties of the top 30 m of the soil column, disregarding the characteristics of the deeper geology. Six site categories have been defined on the basis of average shear wave velocity.

The upper geological layer (soil) has been classified for the whole country using type of bedrock, lithology, and age of the surface layers, as well as some other characteristics taken from detailed geological maps of Yemen Geological Survey on 1 : 200 k scales (Figure 9). The classification follows the NEHRP soil classification scheme (Table 3), which considers seven soil classes and their associated shear wave velocities   ranging from hard rocks (soil index 1.0) to soft soil young water-saturated alluvial deposits (soil index 3.5). In the absence of a Yemen-specific relationship between soil indexes and amplification effects, the widely used NEHRP site amplification procedure [39] has been used in this study. This procedure, together with the RMSI soil classification scheme, provides the amplification factors to scale the ground motion from one soil condition to another. Nonlinear two-dimensional soil amplification factors modified from Choi and Stewart [40], having nonlinear multipliers based on the level of ground motion (PGA) and averaged soil index assigned for a given location, have been used in the study as site-dependent amplification factors.

8. Verification of Ground Motion Estimation Scheme

Before probabilistic hazard calculations, we checked the applicability of the developed strong motion prediction scheme for estimation of ground motion parameters along the territory of Yemen. The intensity data obtained during historical earthquakes (11 January 1941; 13 December 1982, and 13 October 1991, see Figure 10) were used for comparison with the model output. Ground-motion footprints in terms of PGA values under rock site condition were calculated using the averaged values resulting from accepted attenuation equations. The source-to-site distance was estimated as the shortest distance between the fault plane and centroid of a grid/population agglomeration area. Using the (soil index) map, the average values within the grid/population agglomeration area were estimated with GIS (Geographical Information System) overlay analysis. Based on the average soil index value and the estimated PGA values (rock site), two-dimensional soil amplification factors were computed and applied to the rock PGA values. The soil-dependent PGAs were converted into MSK intensities using the PGA-MSK relation shown in Table 4.

The Sa’dah Earthquake of January 11, 1941
This earthquake (epicenter 16.4°N and 43.5°E, intensity MSK VIII) is considered amongst the largest earthquakes in the recent history of Yemen [9]. It caused 1,200 deaths, injured 200 people, and damaged 1,700 houses, out of which 300 were destroyed, 400 were damaged beyond repair, and the rest received minor damage [7]. The magnitude of the earthquake was estimated as 6.5 [2].

The Dhamar Earthquake of December 13, 1982
A large area in Dhamar province of Yemen was rocked by a destructive earthquake (, epicenter 14.7°N and 44.2°E) on 13 December 1982. This moderate shallow event that occurred in a densely populated region about 70 km south of Sana’a resulted in 2,500 deaths and injured 1,500 people [7, 8]. More than 70,000 dwellings were damaged and 500,000 people were affected. Most casualties occurred in the highly populated villages that were having rubble stone masonry and unburned brick houses. Such houses suffered heavy damage in the form of large and deep cracks in walls and the collapse of the outer and inner masonry walls that resulted in partial or complete collapse of the whole construction.

The Al-“Udain and Hazm Al-“Udain (Ibb) Earthquake of November 22, 1991
The November 22, 1991 earthquake ( 4.5, epicenter 13.9°N, 44.1°E) caused widespread damage to housing and infrastructure in Yemen. The most affected areas were the Al-“Udain and Hazm Al-“Udain districts in the Ibb region. Preliminary surveys indicate that 7,150 buildings were damaged and 1,578 were destroyed [41]. However, the reinforced concrete structures in the epicentral area suffered no structural damage indicating the high vulnerability of traditional buildings, especially those constructed on slopes without proper foundation [7].
Obviously, it is not possible to expect a perfect agreement between the available intensity observations and the results of modeling. However, as can be seen from Figure 10, the general features of the regional earthquake effect (epicentral intensity and intensity attenuation) are reproduced quite adequately by the applied scheme of ground-motion modeling.

9. Results of Probabilistic Seismic Hazard Analysis

The following scheme has been applied for the analysis. Ground motion parameter from the stochastic event was calculated at every particular point of the studied territory (centroid of population agglomeration areas, and 10 km grid) using selected ground-motion prediction models and the soil-dependent amplification factor. Every stochastic event is characterized by a return period of occurrence, which depends on the characteristics of seismicity of the given source zone and the magnitude of the event. Thus, the same return period is assigned to the level of ground motion generated by the event. The complete master database of median ground shaking values from all stochastic earthquakes with their rates of occurrence at all locations considering their average site conditions has been prepared. This master database has been used for (a) development of probabilistic seismic hazard maps, that is, return period of particular PGA threshold values at the individual sites, and (b) seismic zoning of the country.

The hazard maps in terms of PGA and MSK intensity with a 39% probability of being exceeded in 50 years (100-year return period) and with a 10% probability of being exceeded in 50 years (475-year return period) are shown in Figures 11 and 12. Note that the site effect has already been incorporated into the hazard assessment. As can be seen, the western and southern parts of Yemen are characterized by the highest level of seismic hazard and the hazard is low across the central and the eastern parts. The largest contributors to Yemen seismic hazard are events from the West Arabian Shield seismic zone.

A seismic zoning scheme for Yemen has been suggested based on hazard values calculated for the 475-year return period. This is the first attempt in developing seismic zoning of the country. The territory of Yemen has been subdivided into 3 seismic zones (Figure 13). By definition, each zone is characterized by a constant hazard, quantified by a different value of the reference peak horizontal ground acceleration as shown in Table 5. Yemen capital area and the Dhamar area are located inside the zone with the highest level of hazard (zone I). Note that the Dhamar area had experienced the destructive earthquake () of 13 December 1982.

Although expected seismic hazard values are calculated at over 30,000 rural and urban population agglomeration areas in Yemen, the hazard values in 30 major cities of Yemen in terms of PGA for 100- and 475-year return periods are presented in Table 6 together with information on the city population from the Yemen 2004 Census.

10. Conclusions

This study is a part of the “The Yemen National Probabilistic Risk Assessment” project funded by The World Bank Global Facility for Disaster Reduction and Recovery (GFDRR). The broad objectives of the project were to apply the principals of probabilistic risk assessment for assessing risks from earthquake, flood, storm surge, tsunami, and landslide in order to develop an analysis of risk exposure and financial response capacity for Yemen. The main output of the project is (i) analysis of the applications of catastrophic risk modeling for hazard risk management in Yemen and (ii) a natural hazard risk atlas for Yemen. The outcome of research described in this paper is a model for probabilistic seismic hazard assessment, which provides necessary inputs for earthquake loss estimation and risk planning on the country level. It is necessary to note that no building code has been developed so far in the country and the constructed seismic zonation scheme may be used as the basis for the code.

We realize, of course, that there are some limitations and shortcomings in our study. We do not consider epistemic uncertainty in all the inputs to the hazard assessment, which usually is incorporated using the logic tree approach, as well as aleatory uncertainty in PGA-intensity conversion. The choice of the attenuation relations used in this study may have some regional biases. The developments of regional ground motion prediction models, as well as the testing of applicability of the recently proposed generalized worldwide ground-motion models, are a matter of high importance. The applied procedure for soil modification requires additional study. For practical calculations, requiring the estimation of the median hazard and loss values, these simplified approaches have been applied.

Acknowledgments

The authors are grateful to the Risk Modeling and Insurance Group, RMSI Pvt. Ltd., Noida (India) for providing necessary facilities. Financial assistance from the World Bank Grant (Contract no. 7149613) is thankfully acknowledged. The authors wish to express their gratitude to Pushpendra Johri. Upamanyu Mahadevan and other members of the Risk Modeling and Insurance Group, RMSI Pvt. Ltd. for their critical reading of the paper. They acknowledge thoughtful comments and suggestions by two anonymous reviewers, which enhanced the quality of paper significantly.