Estimation of VS30 using shallow depth time-averaged shear wave velocity of Rawalpindi–Islamabad, Pakistan

Abstract The time-averaged shear wave velocity of top 30 m (VS30) is the most commonly used parameter to classify a site and evaluate its amplification characteristics for the seismic design. The in-situ seismic tests must be performed up to a depth of 30 m for obtaining the shear wave velocity (VS) profiles to estimate VS30. It is intimated that, in most of the cases, the measured VS profile does not extend up to 30 m due to numerous reasons including limitation of testing techniques and unfavorable field conditions. Since, the measurements of VS30 are unavailable for the majority of Pakistan and the world, the local geology and topographic slope or its combination are used to estimate VS30. However, there is no field-based validation of the estimated VS30 is performed in Islamabad–Rawalpindi region, the proxy-based estimation may lead to unrealistic results. To accommodate this, region specific extrapolation methods are developed. This study develops an empirical data-driven function of VS30 from shallow VS profiles by correlating VS30 with the time-averaged VS to depths less than 30 m. In this regard, 85 VS profiles are used from Rawalpindi-Islamabad region. A comparative analysis of the proposed procedure is carried out with the published methods. It is revealed that VS30 predicted by the proposed function results in close matches with the data measured in the western United States. In addition, the results indicate that the local geology and topographic slope proxies may not be acceptable for usage in the region due to their greater uncertainty. Finally, a procedure for extrapolating the VS profile from available shallow depth measurements up to 30 m is proposed.


Introduction
The shear wave velocity (V S ) of a site represents the dynamic characteristics of a site.It provides the foundation to assess the site amplification factors in quantifying the local site effects due to earthquake shaking using site specific ground response analysis (Kramer 1996).The characteristics of vertically propagated horizontal shear waves from bedrock are significantly influenced by the top 30 m strata (Anderson et al. 1996;Kuo et al. 2011).Therefore, the time-averaged shear wave velocity of the top 30 m (V S30 ) is widely utilized as a proxy to evaluate the variation of the ground motion as it travels through the soil medium to surface (Boore et al. 1997).V S30 is the primary parameter used in various site classification systems (BCP 2007;ICC I 2015).The building code of Pakistan (BCP 2007) is based on the National Earthquake Hazard Reduction Program (NEHRP) provisions (Council BSS 1997) and characterize the sites into five categories listed in Table 1.Furthermore, V S30 is commonly used in the development of site amplification models (Aaqib et al. 2021;Chiou and Youngs 2014;Seyhan and Stewart 2014;Harmon et al. 2019;Tran et al. 2021) which are subsequently employed in the ground motion prediction equations (GMPEs) (Campbell and Bozorgnia 2008) and their ranking and testing for seismic hazard analysis (Bindi et al. 2006;Nizamani and Park 2021) and also for the site-specific response analyses (Nguyen et al. 2020;Tran et al. 2021;Bhusal et al. 2022;Gaha et al. 2022) and microzonation studies (Thompson et al. 2010;Khan and Khan 2018).
A considerable challenge lies in the measurement of V S profiles for regions where the economic conditions are constrained, such as Pakistan (Sadiq et al. 2021).The availability of measured V S profiles up to 30 m is important for characterization of the site to perform a seismic design.However, in most of cases, the profile does not extend up to 30 m due to numerous circumstances including unfavorable field conditions and limitation of the testing equipment (Sun et al. 2005;Sun 2015;Aaqib et al. 2020;2022).In such cases, the V S profiles are empirically derived from the standard penetration test (SPT) data (Manandhar et al. 2018;Adeel et al. 2020;Adeel et al. 2022).Even the availability of the SPT data up to 30 m is not always warranted (Lodi et al. 2015).
Recently, a number of studies have been carried out to estimate V S30 based on various proxy methods such as topography, geology, and geomorphology.Wald and Allen (2007) and Allen and Wald (2009) proposed a procedure to compute V S30 based on the topography of the region in the absence of geological and geotechnical site information.Estimating V S30 as a proxy from horizontal to vertical spectral ratio (HVSR) has also been proposed and widely used (Pokhrel et al. 2019;Yaghmaei-Sabegh and Rupakhety 2020).
Wills and his colleagues (Wills et al. 2000;Wills and Clahan 2006;Wills et al. 2015) (2018) developed seismic site characterization map and classification maps by correlating instrument based V S30 measurements with geological units for northern Pakistan.Given the lack of in situ V S30 observations, these regional-scale proxies-based V S30 calculations can cause large amount of uncertainties.Moreover, due to site characteristics and high risk of earthquake in Islamabad-Rawalpindi region, a more precise V S30 calculation method is needed to evaluate seismic site conditions.
In cases where the measurements do not extend up to 30 m, V S30 has been empirically correlated with shallower time-averaged shear wave velocities and also to extrapolate V S profile up to 30 m. Boore (2004;hereafter referred to as B04) proposed two methods to evaluate V S30 from shallow V S profiles.A total of 135 boreholes from California were utilized for this purpose.The first method assumed a constant velocity from the shallowest depth of available V S profile to 30 m whereas in the second method, the correlation between V S30 and mean V S to the shallower depths was proposed, defined as follows: where V SZ represents the time-averaged V S from surface to the shallow depth (z) whereas a and b are depth dependent regression coefficients.Cadet and Duval (2009) used Eq.(1) proposed by B04 and investigated the correlations between V S30 and average V S to five representative depths of 5, 10, 20, 50 and 100 m.A total of 504 V S profiles from Japan and 22 V S profiles from Europe were utilized.Boore et al. (2011) implemented Eq. (1) to investigate the relationship between V S30 and averaged V S to five representative depths of 5, 10, 15 and 20 m.A total of 638 V S profiles from KiKnet database of Japan, 228 V S profiles from Turkey, 135 V S profiles from California and 21 V S profiles from Europe were utilized.They concluded that the empirical correlation of B04 may not be appropriate for the Japanese sites, which led to the modification of Eq. ( 1) for that region.Sun (2015; hereafter referred to as S15) used 72 V S profiles from Korea to propose a method for calculating V S30 from shallow V S profiles.V S30 was extrapolated from shallow V S profiles by correlating it with the mean V S up to the five representative shallow depths, employing the following function: Furthermore, a function for extrapolating the V S from shallow depths to 30 m was also proposed, defined as follows where z c is the depth at which V S is calculated and V S (z c ) is defined as follows where, z c(cut off) and V s(cut off) indicates the depth of the deepest V s available and V s at that depth, respectively.Mahmood et al. (2016) performed 1D site response analysis of collapsed Margalla Tower in Islamabad during 2005 Kashmir Earthquake.They used SPT measurements of the site and established V S profile using the Lee (1990) correlation.
The site was characterized as Sc according to BCP (2007).A total of four ground motions were used from Pacific Earthquake Engineering Research Center (PEER) database with magnitude ranging 7-7.6 with closest distance.It was concluded that the amplification factors using all the four-input ground motions were amplified up to 2.25, overestimating the design values recommended by BCP (2007).Lodi et al. (2015; hereafter referred to as LEA15) employed a method similar to B04 to compute V S30 from shallow velocity profiles.A total of 140 SPT profiles from the Karachi region of Pakistan with z > 30 m were utilized.However, the V S profiles were estimated from SPT profiles using the empirical correlations.It was reported that the use of SPT profiles to calculate V S adds uncertainty in the conversion (Aaqib et al. 2018).Wang and Wang (2015; hereafter referred to as WW15) used a suite of dataset comprising of 135 boreholes from California and 594 boreholes from KiK-net in Japan to develop and investigate the reliability of proposed method to calculate V S30 from known shallow depth time-averaged shear wave velocities at two depths.A simple extrapolation/interpolation method was proposed to calculate Vs z at a depth.If z is set to 30 m, the equation becomes as follows: where V s(z1) and V s(z2) are time-averaged shear wave velocities at depths z 1 and z 2 , respectively.The limitation of this study is that the whole V S profile is required.Zhou et al. (2021; hereafter referred to as ZEA21) compiled a dataset of 8831 boreholes in China and evaluated several published models based on the dataset to develop a parametric model for estimating V S30 from shallow borehole profile.Similar to WW15, ZEA21 also has the limitation that if the V S profile is not available, the model becomes inapplicable.Midorikawa and Nogi (2015) was one of the evaluated models by ZEA21 and was used to develop the M-Pt model in their study.The functional form used to develop M-Pt model is as follows: where, V Sz30 is defined as follows: where a and b are model specific regression coefficients.
A specific interest of this study is the Rawalpindi-Islamabad region, located in the northern part of Pakistan.The northern part poses a significant threat to the seismic hazard because it lies in the vicinity of the seismically active fault Main Boundary Thrust (MBT).Rapid economic and infrastructure growths has taken place in the capital area of Pakistan, further increasing the societal risk to a devastating earthquake such as the 2005 Kashmir event (moment magnitude ¼ 7.6).This earthquake caused more than 80,000 casualties and 100,000 injuries, while leaving 2.8 million people homeless (Durrani et al. 2005).Some severe repercussions include the collapse of multistory Margalla Towers, Islamabad, located over 80 km from the epicenter.The Building Code of Pakistan (BCP 2007) classifies capital region as Zone 2B, for which the peak ground acceleration (PGA) ranges from 0.19 to 0.26 g.The subsurface geology of the northern part of Pakistan is quite different than southern part of the country (Khan and Khan 2018).As a result, the model proposed by LEA15 cannot be directly used for site characterization except for the host region.
This study proposes a method to predict V S30 from shallow V S profiles in the Rawalpindi-Islamabad region.A total of 85 V S profiles were collected from Rawalpindi-Islamabad. The least squares regression is performed to correlate V S30 with the mean V S at five representative shallow depths.Additionally, a method for extrapolating the Vs profile up to 30 m from shallow V S profiles in the study region is also developed.Furthermore, a comparative study is conducted with the selected published procedures to assess the effectiveness of the proposed method.

Local geology and site profiles
The city of Islamabad is the capital of Pakistan and is located at 33.43 N and 73.04 E, at the brink of Potwar plateau.Rawalpindi is a twin city linked to the capital.The Islamabad region is a section of the sub-Himalayas as reported by Gansser (1964), defined by the Potwar plateau and salt ranges in Pakistan (Gee and Gee 1989).The Rawalpindi-Islamabad metropolitan area occupies one of the chief Himalayan boundary faults (Khan and Khan 2018).Quaternary sediments dominate most of the region.The area consists of the gravelly deposits which are exhibited all along the Margalla hills (Khan et al. 2020).The seismic hazards in the Rawalpindi-Islamabad metropolitan area are substantially connected with various faults in the vicinity of the region, consisting of the Main Boundary Thrust (MBT), the Panjal, Margalla, Hazara, Murree, Jhelum and Mansehra faults (Monalisa et al. 2007;Yong et al. 2008).
In this study, the V S profiles documented in Khan et al. (2020) at 85 locations covering the Rawalpindi-Islamabad region are used.The testing procedures are explained in Khan et al. (2020).The testing locations and site characteristics are shown in Table 2. Figure 1 presents the location of sites where tests were performed.It also shows that the locations of bore holes cover the entire study area.Figure 2 shows V S30 of each soil profile.The sites classified according to BCP (2007) are also presented.A total of 36 profiles falls under site class C, 47 are classified as site class D whereas 2 sites are classified as site class E. V S30 of the profiles is calculated as follows where h i and V Si are the soil layer thickness and shear wave velocity for layer i, respectively.V S30 of the site profiles range from 148 to 554 m/s.Only the profiles with depth up to 30 m are used.Figure 3 shows the distribution of V S profiles among   3 provides a summary of the geology of the Islamabad-Rawalpindi region.This table contains information on geological formation in addition to the mean and median values of V S30 and BCP (2007) site classes.Additional information on the geology of this region may be obtained in the study of Khan et al. (2020).It has been intimated that the measured V S30 can only classify two sites (site classes C and D) within the restricted band of V S30 values.This may be due to the existence of a variety of geological formations as well as a little variation in V S30 values.It is also important to note that the estimation of V S30 cannot be strictly confined to the surficial geology since there is not a distinct transition in the values of V S30 .

Prediction of V S30 using shallow V S profiles
In this section, the correlations of V S30 with V SZ at various shallow depths are derived.In addition, the function development is explained in detail.
The correlations to predict V S30 at shallow depths require two variables: z and the corresponding V SZ .Five representative depths of 5, 10, 15, 20 and 25 m are selected to derive the possible correlations between V SZ (V S05 , V S10 , V S15 , V S20 , V S25 ) and V S30 .These correlations are then generalized and can be used to calculate V S30 from any shallowest depth.Using the collected V S profiles, the V SZ values at the representative depths are calculated.Least squares linear regressions are performed to correlate V SZ with V S30 .Figure 5 illustrates some examples of the correlations between V SZ (V S10 , V S15 , V S20 , V S25 ) and V S30 .It is demonstrated that the coefficient of determination (R 2 ) increases as the depths get closer to 30 m which confirms a stronger correlation of V SZ with V S30 at higher depths.A linear trend is observed which is in line with the studies of B04 and S15.Table 4 summarizes the generic functional form proposed in this study along with the associated R 2 values.It should be noted that the R 2 values computed from the correlations in this study are slightly higher than the published studies, ranging from 0.70 to 0.998.This reveals that proposed functions fits the data well and follow the measurement trend.The functions proposed in this study are similar to S15, as defined in Eq. ( 2) but different than B04 and LEA15.
Figure 6 depicts the correlations between the five selected V SZ and V S30 .At lower V S30 values (V S30 < 200 m/s), the predictions using V SZ at all representative depths are demonstrated to converge.This indicates that the difference between V S30 and V SZ is minimal at soft sites, whereas it increases as the sites become stiffer.In general, the overall trend of these correlations yields a slope that has a propensity to decrease as V sz decreases.These slopes are then analyzed to develop a general function for estimating the V S30 from time-averaged shear wave velocity at any shallower depth.The coefficients of the representative V SZ and V S30 correlations listed in Table 4 are illustrated in Figure 7 at corresponding depths.It is to worth mention here that the slopes of five selected V SZ and V S30 relations are taken as coefficients for average V S (C s ) and then nonlinear regression is performed with respect to z. Figure 7 reveals that Cs increases proportionally with z, C s ¼ 0.67 at z of 5 m and ultimately reaches unity at 30 m. C S and z are strongly correlated and hence provide the basis for developing a generalized empirical correlation between V SZ and V S30 .This trend is in line with the study of S15.The derived correlation is a modified version of the equation proposed in S15 and can be used to calculate V S30 at any shallow depth using the time-averaged shear wave velocity at depth.It is defined as follows: where C S is defined as follows: CS ¼ 0:4643 z 0:2239 (10) Figure 8 compares V S30 estimated by the function proposed in this study with the predictions of V S30 by B04, S15, LEA15, WW15, and ZEA21 at z of 10, 15, 20, and 25 m.As explained in the introduction of this article, WW15 requires time-averaged shear wave velocity at two depths of the borehole, therefore, z 1 was set to 5 m because it was reported to result in a minimum error.For ZEA21, the M-Pt model developed in their study was used in this study because of its better performance, as reported in ZEA21.A comparative analysis revealed that the V S30 values predicted by B04 and ZEA21 results in a close match with the V S30 values estimated in this study.This trend is consistent for all the representative depths considered.S15 results in the highest prediction of V S30 among all the studies presented here.LEA15 correlation over predicts for V S30 < 250 m/s and under predicts otherwise.This demonstrates that the site characteristics of the Rawalpindi-Islamabad region differ from those of Karachi.Hence, the correlation proposed by LEA15 cannot be possibly applied to the region under study.It is to note that all the correlations presented here results in a straight line fit with the exception of the LEA15 model.The empirical model derived by B04 and ZEA21 from a large number of sites are shown to produce surprisingly favorable fits.However, it should be noted that WW15 and ZEA21 requires V S measurements up to 30 m to calculate V S30 from V SZ , which practically reduces the effectiveness of these correlations.Although the B04 and ZEA21 models follow the trend of measured data points, the proposed model exhibits lower residuals as compared to both the models.Figure 9 illustrates the variation of Root Mean Square Error (RMSE) with depth from all the correlations.For each depth z, the V S30 was predicted from Vsz and the RMSE was defined as follows: RMSE reduces for all the correlations as the depths approach 30 m.It is demonstrated that S15 results in the highest RMSE across all depths.RMSE from WW15 and LEA15 is comparable for z ! 10 m.RMSE from the proposed correlations is considerably lower at z ¼ 5 m, whereas, it is comparable with B04 and ZEA21 for Table 4. Correlations between V S30 and V SZ at five representative depths.
z ! 10 m.However, the RMSE produced by this study is lowest among all the studies under investigation across all depths.
The efficacy of the proposed method is assessed by using the ratio of estimated V S30 and measured V S30 (V S30 predicted/V S30 measured).Figure 10 plots the V S30 ratio for all the test locations at the depths of 5, 10, 15, 20, and 25 m.The largest dispersion is observed for the depth of 5 m but the error is not significant.A maximum standard  deviation of 0.12 was calculated among the ratios at all the depths which demonstrates inconsiderable error in the predictions from the proposed correlations.The dispersion decreases as the depth increases.At the depths of 20 and 25 m, most of the points fall on the unity line, demonstrating the efficiency of the proposed method.

Prediction of V S profiles from shallow depths to 30 m
As mentioned earlier, V S profile up to the bedrock or up to 30 m is needed for quantification of ground response under an earthquake loading.In this section, a method to generate the V S profile up to 30 m for the Rawalpindi-Islamabad is developed.A mean V S profile from all the 85 V S profiles used in this study is developed at 1 m intervals up to 30 m.A second order polynomial function similar to that used in S15 is then fitted to the mean V S profile.The following second order polynomial function is proposed to provide V S at any depth up to 30 m, because it was reported that this function provides sufficient fit to the data (Boore et al. 2011): The fitted curve along with the mean V S profile data points is shown in Figure 11.The function is correlated from depths from 5 to 30 m.
To generalize the expression for acquiring V S from any shallow depth for any shallow profile, the expression should describe the final measured depth and V S at that depth.The function is defined as follows: V S (z c ) is defined as follows: Figure 12 illustrates a representative profile measured from field up to 15 m.The function shown in Eq. ( 13) is used to construct the profile from 15 m onwards.V S profile assuming constant velocity from the shallowest available depth up to 30 m is also shown along with the representative measured profile.The proposed expression is demonstrated to effectively extrapolate the profile up to 30 m.The percentage difference between the extrapolated and the measured profile ranges from 2% at 15 m to 9% at 30 m.It is worth mentioning that the use of extrapolated function provides  13) fitted to the mean shear wave velocity profile for depths 5-30 m. much better estimates than assuming the constant V S profile from the point where the field testing is discontinued and up to 30 m.It should be noted that the uncertainties were not accounted for as this was the first attempt for developing the V S profile up to 30 m from shallow profiles in Pakistan.It should also be noted that the developed function is independent of site class and did not account for all possible soil types or soil conditions.These issues should be considered in the future studies.

Conclusions
The evaluation of ground motion characteristics and site classification requires V S30 estimation.In this study, a total of 85 V S profiles measured up to 30 m from Rawalpindi-Islamabad region are utilized to develop the V S30 prediction function from shallow V S profiles.Additionally, a function for extrapolating V S profile from shallow depths to 30 m is also developed.Moreover, the topography and geology of the region are also taken into consideration for a possible correlation with V S30 .It is revealed that topographic slope and geology proxies show a weak with low confidence relationship with measured data.As a result, the region-specific extrapolation method is developed.
The V S30 values of the available profiles are correlated with corresponding timeaveraged shear wave velocity at five representative depths of 5, 10, 15, 20 and 25 m.Least squares linear regression was performed and the functions are demonstrated to have a high coefficient of determination.These correlations yield a generalized function to determine the V S30 from the time-averaged shear wave velocity at shallower depths.A comparative analysis of the proposed procedure is carried out with the published methods.It reveals that the V S30 predicted by the function used in this study results in a close match with those in the western United States and China.
We also analyze the performance of five existing functions and proposed relation with measured data using RMSE.The empirical relation developed in this study outperforms all five models for all depths.The S15 shows the highest RMSE throughout the depths whereas ZEA21 exhibits low errors among other functions.B04 has comparable performances with ZEA21 at all depths except z < 10 m.LEA15 and WW15 show subpar performance.
To extrapolate the V S profile from shallower depths to 30 m, average V S profile of all the sites in this study is used.A quadratic function is used to fit the representative mean V S curve.This function is then utilized to propose a generalize function for extrapolating V S from shallow depths.A comparison with the constant V S extrapolation method reveals that the V S extrapolated using the function proposed in this study results in a realistic V S profile.In general, the empirical function developed for V S30 , and the V S extrapolation method may be used to quantify seismic hazard risk in the region.

Disclosure statement
No potential conflict of interest was reported by the authors.

Figure 1 .
Figure 1.Location points for data acquisition in the Rawalpindi-Islamabad metropolitan area considered in this study.

Figure 2 .
Figure 2. V S30 distribution of the 85 measured profiles as per BCP (2007) site classification system.

Figure 4 .
Figure 4. Correlation of measured Vs 30 with the GMTED2010 digital elevation model (a) topographic slope, and (b) elevation.

Figure 7 .
Figure 7. Correlation between depth (z) less than 30 m and coefficient for average V S (Cs) proposed in Eq. (10).

Figure 8 .
Figure 8.Comparison of the V S30 predictions in this study with the selected published studies in the literature at depths: (a) 10, (b) 15, (c) 20, and (d) 25 m.

Figure 9 .
Figure 9. Variation of RMSE with depth for the correlations under investigation.

Figure 10 .
Figure 10.Variation of predicted and measured V S30 ratio at all the test locations from: (a) V S10 , (b) V S15 , (c) V S20 , and (d) V S25 .

Figure 11 .
Figure 11.Polynomial function proposed in Eq. (13) fitted to the mean shear wave velocity profile for depths 5-30 m.

Figure 12 .
Figure 12.An example of extrapolation of V S profile using the function proposed in Eq. (13).

Table 1 .
Site classification system in Building code of Pakistan (BCP 2007).
Li et al. (2019))al units to estimate V S30 values from the geologic map in California, United States.Xie et al. (2016)employed surface geology and Quaternary depth to the Beijing Plain to develop a V S30 prediction model.Li et al. (2019)took geologic age and grain size into account while constructing a China site categorization map.Recently, Shafique et al.

Table 2 .
Testing locations and characteristics of sites (modified afterKhan et al. 2020).