A Generic Shear Wave Velocity Proﬁling Model for Use in Ground Motion Simulation

: This study presents a generic model for constructing shear-wave velocity (V S ) proﬁles for various conditions that can be used for modeling the upper-crustal modiﬁcation e ﬀ ects in ground motion simulations for seismic hazard analysis. The piecewise P-wave velocity (V P ) proﬁling model is adopted in the ﬁrst place, and the V S proﬁle model is obtained by combining the V P proﬁling model and V S / V P model. The used V S / V P model is constructed from various ﬁeld measurements, experimental data, or CRUST1.0 data collected worldwide. By making the best use of the regionally / locally geological information, including the thickness of sedimentary and crystalline layers and reference V S values at speciﬁc depths, the V S proﬁle can be constructed, and thus the ampliﬁcation behavior of V S for a given earthquake scenario can be predicted. The generic model has been validated by four case studies of di ﬀ erent target regions world around. The constructed proﬁles are found to be in fair agreement with ﬁeld recordings. The frequency-dependent upper-crustal ampliﬁcation factors are provided for use in stochastic ground motion simulations for each respective region. The proposed V S proﬁling model is proposed for region-speciﬁc use and can thus make the ground motion predictions to be partially non-ergodic.


Introduction
Seismic hazard analysis for the regions of low-to-moderate seismicity can be very challenging for the acute scarcity of strong-motion data in these regions, and such regions include southeastern Australia (SEA) continent, eastern North America (ENA), and southeastern China (SEC). A realistic modeling approach for seismic hazard analysis purpose in stable continental regions (SCRs) and other intraplate regions (which are typically considered as low-to-moderate seismicity regions) can be achieved by combining the earthquake source effects, path effects, upper-crustal modification effects, and site conditions through stochastic simulations [1][2][3][4].
The regional path factor (which means path effect model) and local upper-crustal factor (which means the upper-crustal modification effect model and it is different from the widely used site factor, site-transfer function, site effect model, etc.) are of fundamental significance for modeling the ground motions in low-to-moderate seismicity regions with possibly low uncertainties and variabilities [5,6]. In practice, the regional path factor can be identical within a tectonic region in ground motion modeling, and such path factor can be representative of the characteristics of ground motion attenuations within the whole region (with the assumption that the path factor is ergodic within the region). A piece of good evidence for this statement is that most existing ground motion has been adopted nowadays. For the two reasons mentioned above, there is a requirement to propose a more rigorous and comprehensive V S profiling model. The CRUST1.0 database is a global crustal model database at 1 • × 1 • resolution, which is an update of CRUST5.1 and CRUST2.0. The new model database incorporates an updated version of the global sediment thickness. The principal crustal types of the new model database are adopted from CRUST5.1, and the additional crustal types mark specific tectonic settings, such as continental rifts, continental shelves, and oceanic plateaus. In contrast to older models, the function of crustal types in the new model is limited to assigning elastic parameters to layers in the crystalline crust. CRUST1.0 consists of less than 40 crustal types, each of the 1 • × 1 • cells has a unique 8-layer crustal profile where the layers are: (1) water; (2) ice; (3) upper sediments; (4) middle sediments; (5) lower sediments; (6) crystalline upper crust; (7) crystalline middle crust; (8) crystalline lower crust. Parameters like compressional wave velocity (V P ), shear wave velocity (V S ), and crust density (ρ) are given explicitly for each layer (these parameters found the basis of this study). The updated correlation functions between Poisson's ratio and V S , V P and V S , ρ, and V S , are all adopted from Brocher's study which was published in 2005 [20], as this study used the comprehensive dataset to derive the correlation functions and is generally acknowledged by scholars worldwide. More introduction of global crustal models can be found from the paper published by Chandler et al. [3].
The purpose of this study is to propose a more valid and comprehensive V S profiling model, which can be used to construct the region-specific V S profiles for various upper-crustal conditions and can be incorporated into stochastic simulations of ground motions as well as seismic hazard assessment procedures. The model is based on the available geological information which can be obtained from the CRUST1.0 database as well as existing field recordings. The well-known piecewise P-wave model is briefly introduced in Section 2.1. A newly developed V S /V P model as a function of depth (Z) using data from various sources worldwide is put forward in Section 2.2, and 6 frequently encountered cases with the corresponding V S models are identified in Section 2.3. Four case studies from different tectonic regions are conducted in Section 3 to validate the proposed V S profiling model. The final frequency-dependent upper-crustal amplification factor for each region is provided, which can be used in stochastic simulation procedures directly.

Shear-Wave Velocity Modeling
This section introduces the detailed modeling process of the near-source shear-wave velocity (V S ) profile within the upper-crust, which can be used for constructing the frequency-dependent upper-crustal amplification factor. As mentioned above, according to the conclusion obtained by Brocher [20], the estimates of compressional wave velocity (V P ) are reasonably accurate in most cases. However, the estimates of shear-wave velocity (V S ) are not accurate in most cases [21,22]. The reason is that the relationship between V P and V S are always inappropriate. Thus, in this study, V P profiling model would be adopted directly from previous studies in the first place. Then, V S /V P would be modeled using updated recording data collected from various sources. The final Vs profile would then be modeled by combining the V P profile model and V S /V P model.

Compressional Wave Velocity (V P ) Profiling Model
The previous studies show that the compressional wave is closely related to crustal structure, depth, temperature, the geological age of rock formation, and even chemical compositions of the sediments [23][24][25][26]. For seismic hazard analysis purpose, the modern 3D V P profiles obtained from seismic refraction and reflection surveys are not convenient for use in stochastic ground motion procedures, as the purpose of most seismic velocity measurements is for determining site-specific properties, not for developing any specific mathematic expressions [27][28][29]. In this study, all seismic wave profiles will be modeled as the functions of depth for convenient engineering use in stochastic ground motion simulation procedures, which is suggested by other related studies [2,16]. The functional Geosciences 2020, 10, 408 4 of 18 form developed by Chandler et al. [3] is adopted. The piecewise functional form of the V P model is summarized as Equation (1): where Z is the depth in the unit of km; V P0.03 , V PZS , V PZC , V P8 are the reference V P values at the depth of 0.03 km, Z S km, Z C km, and 8 km, respectively; and Z S , Z C are the thickness of the upper sedimentary crustal layer and the total sedimentary crustal layers with the unit of km [3]. The whole function can be illustrated by Figure 1.
where Z is the depth in the unit of km; VP0.03, VPZS, VPZC, VP8 are the reference VP values at the depth of 0.03 km, ZS km, ZC km, and 8 km, respectively; and ZS, ZC are the thickness of the upper sedimentary crustal layer and the total sedimentary crustal layers with the unit of km [3]. The whole function can be illustrated by Figure 1. The validation of this VP model can be found in Chandler et al. [3], and the results showed that the VP model is valid in most cases. Therefore, this VP model is adopted directly.

Modeling of VS/VP
The VS-VP relationship is essential for various engineering applications, including the analysis of reservoir geomechanical properties for studying seismic ground motions. However, VS records are often not available for target regions due to technological limitations. Alternatively, VS can be obtained from empirical prediction equations based on the records of VP. In the BJ97 [2] model, Boore and Joyner obtained the VS profiles generally from two generic approaches: one is the VS data from boreholes in the upper 4 km; another one is through the data of the VP profile and the VS-VP relationship (VS/VP = 1/√3) for the depth larger than 4 km. In this study, the VS profiles will be derived from VP completely for all depths (as a generic approach). Thus, the VS-VP relationship is of extreme importance to obtain a relatively accurate VS profile (assume VP profile is valid). Recording data from multiple sources around the world are gathered for modeling VS/VP. These data are selected randomly world around to make the VS/VP to be as "generic" as possible. Detailed information about data collected for modeling VS/VP is summarized in Table 1. The validation of this V P model can be found in Chandler et al. [3], and the results showed that the V P model is valid in most cases. Therefore, this V P model is adopted directly.

Modeling of V S /V P
The V S -V P relationship is essential for various engineering applications, including the analysis of reservoir geomechanical properties for studying seismic ground motions. However, V S records are often not available for target regions due to technological limitations. Alternatively, V S can be obtained from empirical prediction equations based on the records of V P . In the BJ97 [2] model, Boore and Joyner obtained the V S profiles generally from two generic approaches: one is the V S data from boreholes in the upper 4 km; another one is through the data of the V P profile and the V S -V P relationship (V S /V P = 1/

√
3) for the depth larger than 4 km. In this study, the V S profiles will be derived from V P completely for all depths (as a generic approach). Thus, the V S -V P relationship is of extreme importance to obtain a relatively accurate V S profile (assume V P profile is valid). Recording data from multiple sources around the world are gathered for modeling V S /V P . These data are selected randomly world around to make the V S /V P to be as "generic" as possible. Detailed information about data collected for modeling V S /V P is summarized in Table 1. Thorough regression analysis for the ratio (V S /V P ) is conducted using this dataset. As the data set is collected world around, large uncertainties should be prudently considered during the modeling process. Chandler et al. [3] provided a power functional format for the V S /V P ratio as it is convenient to combine it with the V P model. In this study, the authors will adopt the power functional format again but in a more rigorous way. The detailed modeling procedures are stated below.
Firstly, as Poisson's ratio (denoted as "σ" in this study) is closely related to the V S /V P ratio, it is important to clarify it before the determination of the V S /V P ratio. For shallow depth, σ value depends largely on the depth and saturation of the medium. However, it is suggested by Boore and Joyner [2] that σ can be fixed at 0.25 for the depth larger than 4 km, as the lithostatic pressure at these depths (about 1 kbar) can make most cracks be closed. From the study of Brocher [20], the functional form is expressed in Equation (2), which is the same one shown in Chandler et al. [3], which is obtained from the elastic wave propagation theory in a homogeneous medium: where σ is Poisson's ratio. For the depth larger than 4 km, V S /V P = 0.577 (being 1/ √ 3), which is suggested by the BJ97 model and CLT05 model. According to this study, for the depths between 2 and 4 km, the V S /V P value also keeps constant (shown in Figure 2). To be more explicit, the authors found that below the depth of 2 km, V S /V P = 0.577 can obtain similar regression goodness (with no apparent change in mean residual, which is defined by Equation (3), for the depths larger than 2 km). Therefore, for the depth larger than 2 km, V S /V P is set to be a fixed value 0.577 in this study (Zone C in Figure 2): where δ is the mean residual of all depth (the upper-limit of depth is 50 km if no specific note is given); N is the total number of recording data collected from various sources, which equals to 193 in this study; y i is the field recorded V S /V P value;ŷ i is the predicted V S /V P value, which is obtained from the proposed model in this study.
Different transition depth values (0.1, 0.2 , 0.4, and 1 km) are evaluated. The first segment would only consider modeling the dataset with the depth smaller than the transition depth. The authors found that the mean residual can achieve the lowest value for the data with depth less than 0.2 km. Thus, the transition depth is set to be 0.2 km in this study (shown as Zone A in Figure 2). Lastly, another segment (from 0.2 to 2.0 km) is modeled to connect the two segments modeled in the previous steps, and the result is shown as Zone B in Figure 2. As shown in Figure 2, the proposed tri-linear VS/VP model can represent the collected data (listed in Table 1) with a reasonable level of accuracy (with the mean residual of 0.089). The point should be noted is that this generic VS/VP model is proposed (as a sole function of depth) for convenient use in engineering seismology and other related studies, other parameters like material type, temperature, and pressure are not considered.
The overall VS/VP model is expressed in Equation (4): To evaluate the rationality of the proposed VS/VP model in this study, four existing VS-VP relationship equations obtained from various data recordings globally are adopted in this study for comparison purpose (the equations can be found in Table 2). Note that all the selected equations, Cea85 [42], Han86 [43], B05 [20], CLT05 [3] are used as generic models for global conditions. To illustrate the VS/VP ratio with the change of depth, the generic VS profile of NEHRP (National Earthquake Hazard Reduction Program) B/C site (VS30 = 0.76 km/s) is adopted and VS/VP is computed Secondly, for the depth smaller than 2 km, the general regression technique is performed. Single-segment line (log(V S /V P ) vs. log(Z)) is adopted in the first time, the mean residual value (0.21) indicates the regression can be improved to some extent. In this case, two segments are considered. Different transition depth values (0.1, 0.2, 0.4, and 1 km) are evaluated. The first segment would only consider modeling the dataset with the depth smaller than the transition depth. The authors found that the mean residual can achieve the lowest value for the data with depth less than 0.2 km. Thus, the transition depth is set to be 0.2 km in this study (shown as Zone A in Figure 2).
Lastly, another segment (from 0.2 to 2.0 km) is modeled to connect the two segments modeled in the previous steps, and the result is shown as Zone B in Figure 2.
As shown in Figure 2, the proposed tri-linear V S /V P model can represent the collected data (listed in Table 1) with a reasonable level of accuracy (with the mean residual of 0.089). The point should be noted is that this generic V S /V P model is proposed (as a sole function of depth) for convenient use in engineering seismology and other related studies, other parameters like material type, temperature, and pressure are not considered.
The overall V S /V P model is expressed in Equation (4): To evaluate the rationality of the proposed V S /V P model in this study, four existing V S -V P relationship equations obtained from various data recordings globally are adopted in this study for comparison purpose (the equations can be found in Table 2). Note that all the selected equations, Cea85 [42], Han86 [43], B05 [20], CLT05 [3] are used as generic models for global conditions. To illustrate the V S /V P ratio with the change of depth, the generic V S profile of NEHRP (National Earthquake Hazard Reduction Program) B/C site (V S30 = 0.76 km/s) is adopted and V S /V P is computed based on this profile. The comparison details are shown in Figure 3 for depth ranging from 0.001 to 100 km for all models. Table 2. V S -V P relationship equations used for comparison purpose.

Equations
Citations Geosciences 2020, 10, x FOR PEER REVIEW 7 of 19 based on this profile. The comparison details are shown in Figure 3 for depth ranging from 0.001 to 100 km for all models. Table 2. VS-VP relationship equations used for comparison purpose. As shown in Figure 3, except for the CLT05 model, other models can obtain similar estimates of VS/VP for different depths. It can be seen that the proposed VS/VP by this study matches well with the B05 model, which was developed from comprehensive datasets [20] and widely acknowledged globally. The CLT05 model gives smaller estimates of VS/VP at the depth 0.01 to 4 km, compared with other models. The Cea85 model gives similar estimates of VS/VP at the depth smaller than 2.0 km, compared with the B05 model and the proposed model in this study. The Han86 model gives a similar estimate of VS/VP at the depth of 0.2 km, and larger estimates of VS/VP at other depths, compared with the B05 model and the proposed model in this study. The mean residuals (using Equation (3)) between the proposed model and Cea85, Han86, B05, and CLT05 models are 0.039, 0.047, 0.012, 0.053 respectively, which indicates that the proposed VS/VP model is close to the value given by the B05 model for all depths.

Shear-Wave Velocity Profile Modeling
Based on the results obtained from Sections 2.1 and 2.2, the VS profile can be modeled by combining the VP profile and VS/VP model. The detailed modeling processes are stated below.
For the Zone I and Zone A (Zone IA): For the Zone I and Zone B (Zone IB):  As shown in Figure 3, except for the CLT05 model, other models can obtain similar estimates of V S /V P for different depths. It can be seen that the proposed V S /V P by this study matches well with the B05 model, which was developed from comprehensive datasets [20] and widely acknowledged globally. The CLT05 model gives smaller estimates of V S /V P at the depth 0.01 to 4 km, compared with other models. The Cea85 model gives similar estimates of V S /V P at the depth smaller than 2.0 km, compared with the B05 model and the proposed model in this study. The Han86 model gives a similar estimate of V S /V P at the depth of 0.2 km, and larger estimates of V S /V P at other depths, compared with the B05 model and the proposed model in this study. The mean residuals (using Equation (3)) between the proposed model and Cea85, Han86, B05, and CLT05 models are 0.039, 0.047, 0.012, 0.053 respectively, which indicates that the proposed V S /V P model is close to the value given by the B05 model for all depths.

Shear-Wave Velocity Profile Modeling
Based on the results obtained from Sections 2.1 and 2.2, the V S profile can be modeled by combining the V P profile and V S /V P model. The detailed modeling processes are stated below.
For the Zone I and Zone A (Zone IA): Geosciences 2020, 10, 408 8 of 18 For the Zone I and Zone B (Zone IB): For the Zone I and Zone C (Zone IC): Similarly, the V S model for the Zone IIIA, IIIB, and IIIC are expressed as Equations (8)-(10), respectively: For the transition zone (II): in which: According to the relative positions of Z S and Z C to the identified marker depth 0.2 and 2.0 km, there would be 6 (which is 3!) different cases, which are illustrated in Figure 4.
For the Zone I and Zone C (Zone IC): Similarly, the VS model for the Zone IIIA, IIIB, and IIIC are expressed as Equations (8)-(10), respectively: For the transition zone (II): in which: According to the relative positions of ZS and ZC to the identified marker depth 0.2 and 2.0 km, there would be 6 (which is 3!) different cases, which are illustrated in Figure 4. The detailed functional expressions for the 6 identified cases and their corresponding zones are summarized in Table 3. The detailed functional expressions for the 6 identified cases and their corresponding zones are summarized in Table 3. Table 3. The functional forms of V S profile model in this study. In all cases, Zone II is a transition zone that is used to connect the upper zone and the lower zone. Thus, the functional form is the same, except for the value of exponent "n". In all cases, n is expressed by Equation (12). As the values of V SZC and V SZS are changing for different cases, the value of n will be differing accordingly case by case.

Case Depth Range (km) V S (km/s) Zone
The estimation of Z C and Z S values are essential for constructing reasonable local V S profiles. The detailed three approaches for estimating Z C and Z S values are given below.
(i) If sufficient and reliable field measurements of the crustal V S profile are available, the values of Z S and Z C can be adjusted to achieve a good match with the local field measurements; (ii) If no field measurements can be obtained to determine Z S and Z C , then Z S and Z C can be taken as the averaged thickness of the upper sediment layer and averaged total sediment layers (provided by CRUST1.0), respectively.
To make the overall process of V S profile construction more convenient for use, a MATLAB-based computer program is provided alongside this article. The flowchart of the program is shown in Figure 5, and the program can be downloaded from the link: https://github.com/Y-Tang99/GMSS.  Table 3.

Validation of the Proposed VS Profiling Model
In this section, the VS profiling model proposed in Section 2 will be applied to four different tectonic regions. Local VS profiles obtained from technical measurements at each region will be used for comparison and validation purpose. The selected regions are the Melbourne Region (located in southeastern Australia), St. Louis Metro Region (located in central eastern North America), Hong Kong Region (located in southeastern China), and northern Switzerland (located in central western Europe). These regions are selected from typical intraplate regions with low-to-moderate seismicity. The detailed location information of sampling points used for collecting geological information from CRUST1.0 is listed in Table 4 and shown in Figure 6. These sampling points are selected because the local field measurements of VS profiles are collected from these points. The parameters used for constructing VS profiles are summarized in Table 5.
Melbourne Region locates at Victoria state in southeastern Australia. The field measurements of the 5 selected sites in different suburbs around Melbourne metropolitan area are obtained from the passive seismic investigation technique called the spatial auto-correlation (SPAC) method [44]. A series of surveys were carried out to get the VS profiles down to the depth of around 0.1 km into Silurian mudstone. The St. Louis Metro Region is located at Missouri urban area in central eastern North America (CENA). The filed VS profiles of the 6 sites are determined from the surface wave measurements using the multi-channel analysis of surface waves (MASW) geological technique to the depth of around 0.03 km of upper sedimentary bedrock or surficial material [45]. Hong Kong Region locates at southeastern China. The field VS recordings are obtained from several sources: for the depth up to around 0.1 km, the VS data are collected from extensive borehole data; for the depth up to 1.5 km, the VS data are obtained from the SPAC method and short-period group velocity (T = 0.4 − 1.3 s) dispersion of Rg waves generated by quarry blasts [46]. Northern Switzerland is in central western Europe. The field VS profiles are obtained from array processing of ambient noises [47] as well as from interpretation of seismic refraction and reflection studies [48,49] for the depth up to around 3.0 km. For all selected regions, the VS data for the depth exceeding 4 km are obtained from the web-based global crust model CRUST1.0 [50]. The average field measurement for each region is obtained from the mean value of the field measurements at each depth.
For each selected region, the CLT05 model and BJ model will be adopted for comparison purposes. The overall results are shown in Figure 7, and the residuals (defined as the difference  Table 3.

Validation of the Proposed V S Profiling Model
In this section, the V S profiling model proposed in Section 2 will be applied to four different tectonic regions. Local V S profiles obtained from technical measurements at each region will be used for comparison and validation purpose. The selected regions are the Melbourne Region (located in southeastern Australia), St. Louis Metro Region (located in central eastern North America), Hong Kong Region (located in southeastern China), and northern Switzerland (located in central western Europe). These regions are selected from typical intraplate regions with low-to-moderate seismicity. The detailed location information of sampling points used for collecting geological information from CRUST1.0 is listed in Table 4 and shown in Figure 6. These sampling points are selected because the local field measurements of V S profiles are collected from these points. The parameters used for constructing V S profiles are summarized in Table 5.   To use the generic VS profiling model proposed by this study, the ZS and ZC are required to be determined in the first place. The followed step is to determine which case model (listed in Table 3) is suitable for use in the target region. The next step is to determine the reference VS values at certain depths and the exponential value (n). The final step is to construct the VS profile for the target region. To make it clearer, Melbourne Region is taken as an example: Step (i): ZS and ZC values are determined from the average value of the 5 sampling points (for each sampling point, the detailed ZS and ZC value can be collected from CRUST1.0 database). In this study, ZS = 0.05 km and ZC = 4.0 km; Step (ii): Because ZS < 0.2 < 2 ≤ ZC, the Case 4 VS profiling model should be used in this study (listed in Table 3); Step (iii): Ideally, the reference VS values should be determined from local field measurements, while the VS values collected from CRUST1.0 database can be used if valid field measurements are not available. For the Case 4 profiling model, four reference VS values need to be determined: VSZS, VSZC, VSZI, and VS8. VSZS = 1.33 km/s and VSZI = 1.1 km/s (in this case VSZI = VS0.03), and they are determined from average field measurement. VSZC = 3.3 km/s and VS8 = 3.5 km/s, which are determined from the CRUST1.0 database. "n" is finally determined using Equation (12), which is 0.21 in this study.
Step (iv): Construct VS profile for all depths and eliminate any possible abnormal points by a double check procedure. The final VS profile constructed for the Melbourne Region is shown as the red line in Figure 7a. Figure 7a indicates that all model estimates can get the alignment with the field measurements with a reasonable level of accuracy. The average residual of for the three models (this study, CLT05 model and BJ model) is 0.103, 0.107, and 0.075 respectively, which is shown in Figure 8a.
For the St. Louis Metro Region, the situation is slightly different from other selected regions as this region is located at St. Louis Metropolitan area and the sampling sites are quite close to each  Melbourne Region locates at Victoria state in southeastern Australia. The field measurements of the 5 selected sites in different suburbs around Melbourne metropolitan area are obtained from the passive seismic investigation technique called the spatial auto-correlation (SPAC) method [44]. A series of surveys were carried out to get the V S profiles down to the depth of around 0.1 km into Silurian mudstone. The St. Louis Metro Region is located at Missouri urban area in central eastern North America (CENA). The filed V S profiles of the 6 sites are determined from the surface wave measurements using the multi-channel analysis of surface waves (MASW) geological technique to the depth of around 0.03 km of upper sedimentary bedrock or surficial material [45]. Hong Kong Region locates at southeastern China. The field V S recordings are obtained from several sources: for the depth up to around 0.1 km, the V S data are collected from extensive borehole data; for the depth up to 1.5 km, the V S data are obtained from the SPAC method and short-period group velocity (T = 0.4-1.3 s) dispersion of Rg waves generated by quarry blasts [46]. Northern Switzerland is in central western Europe. The field V S profiles are obtained from array processing of ambient noises [47] as well as from interpretation of seismic refraction and reflection studies [48,49] for the depth up to around 3.0 km. For all selected regions, the V S data for the depth exceeding 4 km are obtained from the web-based global crust model CRUST1.0 [50]. The average field measurement for each region is obtained from the mean value of the field measurements at each depth.
For each selected region, the CLT05 model and BJ model will be adopted for comparison purposes. The overall results are shown in Figure 7, and the residuals (defined as the difference between average field measurements and the model predictions) are shown in Figure 8 for each region.
To use the generic V S profiling model proposed by this study, the Z S and Z C are required to be determined in the first place. The followed step is to determine which case model (listed in Table 3) is suitable for use in the target region. The next step is to determine the reference V S values at certain depths and the exponential value (n). The final step is to construct the V S profile for the target region. To make it clearer, Melbourne Region is taken as an example: Step (i): Z S and Z C values are determined from the average value of the 5 sampling points (for each sampling point, the detailed Z S and Z C value can be collected from CRUST1.0 database). In this study, Z S = 0.05 km and Z C = 4.0 km; Step (ii): Because Z S < 0.2 < 2 ≤ Z C , the Case 4 V S profiling model should be used in this study (listed in Table 3); Step (iii): Ideally, the reference V S values should be determined from local field measurements, while the V S values collected from CRUST1.0 database can be used if valid field measurements are not available. For the Case 4 profiling model, four reference V S values need to be determined: V SZS , V SZC , V SZI , and V S8 . V SZS = 1.33 km/s and V SZI = 1.1 km/s (in this case V SZI = V S0.03 ), and they are determined from average field measurement. V SZC = 3.3 km/s and V S8 = 3.5 km/s, which are determined from the CRUST1.0 database. "n" is finally determined using Equation (12), which is 0.21 in this study.
Step (iv): Construct V S profile for all depths and eliminate any possible abnormal points by a double check procedure. The final V S profile constructed for the Melbourne Region is shown as the red line in Figure 7a. other. This site is characterized as the NEHRP rock site (0.62 km/s < VS30 < 1.5 km/s), which is the same as the Melbourne Region and northern Switzerland, while Hong Kong Region is a typical hard rock site (VS30 > 1.5 km/s) [15]. Repeating Step (i) to Step (iv), and the Case 2 profiling model is suitable for this region. Figures 7b and 8b indicate the model proposed by this study performs slightly better than another two models (CLT05 model and BJ model).
For the Hong Kong Region, the average depth of the sedimentary layer indicates that CRUST1.0 is very small. Considering the quite small VS gradient at shallow depths indicated by the field measurements, the exponential value (which indicates the gradient of VS) in the VS profiling model needs to be small accordingly. However, the exponential value for shallow depths (<0.2 km) is fixed at 0.3297 in all case models listed in Table 5, which is too large to represent the actual VS gradient. Based on these considerations, the ZS value is adjusted to be 0.001 km to make the shallow depth zone located in Zone II (shown in Figure 4). In this case, a variable exponential value ("n") can be obtained and can thus be adjusted to mimic the VS gradient of VS obtained from field measurements. As ZC is determined as 1.0 km from CRUST1.0 database, the Case 5 profiling model should be adopted in this region (ZS < 0.2 < ZC ≤ 2). The results are shown in Figures 7c and 8c, and they indicate that the model proposed by this study performs better than CLT05 model and BJ model, especially at shallow depths (<0.01 m).
Similar to Hong Kong Region, the sedimentary layer depth of northern Switzerland is very small indicated by CRUST1.0, and the field VS gradient is small at shallow depths. The value of ZS is set to be 0.001 km, and the Case 5 profiling model is adopted. Final VS profiles and the corresponding residuals are shown in Figures 7d and 8d, respectively. The results also indicate the model proposed by this study performs better than the CLT05 model and BJ model for this region.  The upper-crust amplification can be obtained from the VS profile and density profile. The amplification factor obtained from the quarter wavelength approximation (QWA) method for each region is provided in this study. The frequency-dependent amplification is computed using Equation (13): where ρ and β are the density and VS near the seismic source; ρ and β are the time-weighted average density and VS within the crust. In this study, the units of ρ and β are g/cm 3 and km/s, respectively. The density profiles (ρ) for each region are determined using the equations provided in  Figure 7a indicates that all model estimates can get the alignment with the field measurements with a reasonable level of accuracy. The average residual of for the three models (this study, CLT05 model and BJ model) is 0.103, 0.107, and 0.075 respectively, which is shown in Figure 8a.
For the St. Louis Metro Region, the situation is slightly different from other selected regions as this region is located at St. Louis Metropolitan area and the sampling sites are quite close to each other. This site is characterized as the NEHRP rock site (0.62 km/s < V S30 < 1.5 km/s), which is the same as the Melbourne Region and northern Switzerland, while Hong Kong Region is a typical hard rock site (V S30 > 1.5 km/s) [15]. Repeating Step (i) to Step (iv), and the Case 2 profiling model is suitable for this regionFigures 7b and 8b indicate the model proposed by this study performs slightly better than another two models (CLT05 model and BJ model).
For the Hong Kong Region, the average depth of the sedimentary layer indicates that CRUST1.0 is very small. Considering the quite small V S gradient at shallow depths indicated by the field measurements, the exponential value (which indicates the gradient of V S ) in the V S profiling model needs to be small accordingly. However, the exponential value for shallow depths (<0.2 km) is fixed at 0.3297 in all case models listed in Table 5, which is too large to represent the actual V S gradient. Based on these considerations, the Z S value is adjusted to be 0.001 km to make the shallow depth zone located in Zone II (shown in Figure 4). In this case, a variable exponential value ("n") can be obtained and can thus be adjusted to mimic the V S gradient of V S obtained from field measurements. As Z C is determined as 1.0 km from CRUST1.0 database, the Case 5 profiling model should be adopted in this region (Z S < 0.2 < Z C ≤ 2). The results are shown in Figures 7c and 8c, and they indicate that the model proposed by this study performs better than CLT05 model and BJ model, especially at shallow depths (<0.01 m).
Similar to Hong Kong Region, the sedimentary layer depth of northern Switzerland is very small indicated by CRUST1.0, and the field V S gradient is small at shallow depths. The value of Z S is set to be 0.001 km, and the Case 5 profiling model is adopted. Final V S profiles and the corresponding residuals are shown in Figures 7d and 8d, respectively. The results also indicate the model proposed by this study performs better than the CLT05 model and BJ model for this region.
The upper-crust amplification can be obtained from the V S profile and density profile. The amplification factor obtained from the quarter wavelength approximation (QWA) method for each region is provided in this study. The frequency-dependent amplification is computed using Equation (13): where ρ 0 and β 0 are the density and V S near the seismic source; ρ Z and β Z are the time-weighted average density and V S within the crust. In this study, the units of ρ 0 and β 0 are g/cm 3 and km/s, respectively. The density profiles (ρ) for each region are determined using the equations provided in [15]. More detailed information about the QWA method can be found in [2]. The amplification for each target region can be found in Figure 9.
Geosciences 2020, 10, x FOR PEER REVIEW 15 of 19 [15]. More detailed information about the QWA method can be found in [2]. The amplification for each target region can be found in Figure 9.

Discussion
Rather than providing an approach of estimating the site indicator VS30, this study intends to provide an approach to obtain the shear-wave velocity profile for modeling the upper-crustal modification effects. There are generally three different methods to estimate the upper-crustal factor [51]. The most widely used method is the reference-rock-site method, which is based on the comparison of recorded motions on the local site to those at an identified reference rock site. The reference rock site is often set to be identical for the whole target region (e.g., WNA). The second method is "H/V" (the ratio of horizontal component motions to vertical component motions) method, which is based on the assumption that vertical component ground motions are not influenced by the local site conditions significantly. Thus, the H/V ratio can be a good indicator of local site condition influence on the horizontal component ground motions [52]. This method is similar with the receiverfunction method initially used for studying the earth upper mantle and crust from teleseismic recordings [53], and many existing GMPEs derive the site factor using this method [5,10,13,54]. However, this method often does not consider the effect of upper-crust (the depth is too small). Moreover, a criticism of the "H/V" method summarized by Atkinson and Boore [5] is that the "H/V"

Discussion
Rather than providing an approach of estimating the site indicator V S30 , this study intends to provide an approach to obtain the shear-wave velocity profile for modeling the upper-crustal modification effects. There are generally three different methods to estimate the upper-crustal factor [51]. The most widely used method is the reference-rock-site method, which is based on the comparison of recorded motions on the local site to those at an identified reference rock site. The reference rock site is often set to be identical for the whole target region (e.g., WNA). The second method is "H/V" (the ratio of horizontal component motions to vertical component motions) method, which is based on the assumption that vertical component ground motions are not influenced by the local site conditions significantly. Thus, the H/V ratio can be a good indicator of local site condition influence on the horizontal component ground motions [52]. This method is similar with the receiver-function method initially used for studying the earth upper mantle and crust from teleseismic recordings [53], and many existing GMPEs derive the site factor using this method [5,10,13,54]. However, this method often does not consider the effect of upper-crust (the depth is too small). Moreover, a criticism of the "H/V" method summarized by Atkinson and Boore [5] is that the "H/V" estimate is largely a measure of Rayleigh wave ellipticity and the "H/V" ratio would be largely controlled by site response when it applied to body waves which are measured from earthquakes [52]. Many existing GMPEs are developed by combining the reference rock site method and "H/V" method [12,13]. The third method is the theoretical wave-propagation modeling of upper-crust modification effect (theoretically modeling method). In this method, two separate effects are considered to model the upper-crustal modification effects: upper-crustal amplification effect of the seismic waves when they cross the boundary between different mediums and the upper-crustal attenuation effect which is associated with the transmission quality of the upper-crustal layers (the attenuation effect is beyond the scope of this study). The two effects mentioned above have been observed from the instrumental records from deep drill-holes in active seismic areas [55]. Chen [51] proved that the third method could be comparable to "H/V" method for various upper-crustal conditions.
In the theoretically modeling method, V S and density (ρ) profiles are essentially required to model the upper-crustal amplification effect. The quarter wavelength approximation (QWA) and the square root impedance (SRI) method [2,56] are used to compute the frequency-dependent upper-crustal amplification factor. This method will be adopted in this study to compute the upper-crustal amplification factor.
Seismic wave velocity (including compressional velocity (V P ) and shear-wave velocity (V S )) model can be obtained from non-invasive seismic techniques, including refraction and reflection surveys [57,58], and invasive seismic techniques, including borehole methods [59]. However, most of the seismic wave velocity structures (or profiles) directly obtained by these approaches are often used for deterministically analyzing seismic hazard purpose [60][61][62]. The so obtained seismic wave velocity structures are not convenient (and often not available in low-to-moderate seismicity regions) for use in stochastic ground motion simulations because no parameterized functional form is given. Therefore, seismic waven profile modeling still requires more attention for convenient engineering application, especially for the regions where stochastic simulations are required for seismic hazard analysis.
With the rapid development of the non-ergodic PSHA, the requirement of the site-specific GMPE modeling becomes more acute. This study mainly provides an approach to minimize the uncertainties of upper-crustal modification effects in the modeling procedures, which would help to reduce the aleatory uncertainty in PSHA and make the hazard curve more reliable. Figures 7-9 indicate that, the shallower part of V S profiles critically influences the estimates of upper-crust modification effects as they contain large uncertainties, which means the ground motions at higher frequencies would be more impacted. In practical use, users should pay more attention to the shallower part of the constructed V S profile for the target region and make sure the uncertainties are fully considered in ground motion simulations, and in any following procedures such as PSHA or dynamic structural response analysis using the simulated accelerograms.

Conclusions
This study proposed a generic V S profiling model which can be used for stochastic simulation of ground motions and seismic hazard assessment purposes. The proposed model is mainly designed for rock site conditions and should be avoided for any arbitrary applications. A comprehensive V S /V P relationship, as a function of depth, has been proposed in the first step based on various field recording data collected worldwide. The region-specific V S profile could be constructed using the functional formats provided in Table 3 and the algorithm shown in Figure 5. The geological information needed for V S profile construction includes the thickness of the upper sedimentary layer and the total thickness of all sedimentary layers (collected from CRUST1.0 with the longitude and latitude of the sampling points in the target region); and the reference V S values at certain depths (collected from field recordings or CRUST1.0). The detailed steps of constructing V S profiles have been provided. The validation of the proposed model has been made by 4 case studies in different tectonic regions in the world.