Determining 1D velocity model from local earthquake data in the South Hangay region, central Mongolia

One dimensional (1D) velocity models are still widely used for computing earthquake locations at seismological centers. The location accuracy of an earthquake strongly depends on the velocity model used to compute the location. In the past, the local velocity model developed for the Hangay region was lacking precision due to insufficient data. Within the framework of the “Intracontinental Deformation and Surface UpliftGeodynamic Evolution of the Hangay Dome, Mongolia, Central Asia” project [15], 72 seismic Broadband stations network were deployed in the Hangay Dome. This gives us an opportunity to estimate the crustal velocity structure of the South Hangay region using recorded local earthquake data. For this purpose, available velocity models for the South Hangay region have been re-evaluated. By simultaneous invertion Pand S-wave arrival times using VELEST algorithm, we estimated minimum 1D velocity models, station corrections, hypocentre locations, and origin times for the south Hangay region. Consequently, 1D crustal velocity model is proposed for the South Hangay region. This new model is expected to improve the accuracy of the routine hypocenter determination and as initial reference models for seismic tomography study.


INTRODUCTION
South Hangay region, the study area is located in the southern part of the Hangay Mountains in central Mongolia. The Hangay Mountains cover an area of more than 200,000 km 2 of uplift in central Mongolia with a maximum altitude of 4021 m at the summit of Otgon Tenger Uul. It has a dome shape; the mountain is also called "Hangay Dome". The Hangay Dome occurrs in a kinematic transition between predominantly compressional deformation to the south and extensional deformation to the north. _______________________________________________________________________________

PMAS
The northern side of the plateau and adjacent parts of the Siberian Craton are dominated by the active Baikal and Khuvsgul intracontinental rifts [1,2,3]. The Hangay Dome is bound by three major strike-slip fault systems. These include the Bulnai, the Gobi-Altai and the Mongolian Altai fault system to the north, south, and west, respectively [3,6]. The Gobi-Altai and Bulnai systems are dominated by E-W left-lateral strike-slip displacement faults, which are seismically active. Mongolian Altai fault system is dominated by right-lateral strike-slip faults ( Figure 1A). Together these fault systems accommodate the majority of north-south shortening in western Mongolia and counterclockwise vertical axis rotation of crustal blocks within the Altai Mountains [1,2,10].
The Hangay Mountains are sandwiched between the Bulnai region to the north and the Gobi-Altai region to the south. Few earthquakes have been recorded from the center and southern flank of the dome which is the focus of this study. This study region covered within 45.7°-47.5°N, 97°-101°E coordinates that includes the South Hangay fault system. It extends for ~350 km along the southern slopes of the Hangay [3,4,5] and connects clear evidence of left-lateral-strike-slip faulting in the late Quaternary into a single structure, clustered around the town of Bayanbulag [6]. Farther east, the NW-SE trending Bayankhongor fault, has been described both as a thrust and a normal fault ( Figure 1A and 1B) [1,3].   [5] Proceedings of the Mongolian Academy of Sciences PMAS Due to the lack of recorded earthquake information within the Hangay region in the 20 th century, the area was assumed to be aseismic [2,3,4,5,6]. Except for several areas with high seismic activity during the last century, generally seismic rate is low. Seismicity within the Hangay Dome is particularly sparse [8,10]. However, the seismic activity observed during last century does not necessarily reveal regions with the potential to rupture in large earthquakes. For example, the Hangay region experienced large events 300 to 500 years ago as revealed by the Egiin Davaa normal fault morphology. It is localized in the middle of the dome and could have ruptured in second half of the 16 th century [7] in an earthquake with a magnitude of more than 7 on the Richter scale.
One of the main segments of the South Hangay fault system was activated by a moderate size earthquake with local magnitude Ml=5.4 on 10 March of 2012. This is the most recent strongest earthquake that had occurred in the region. The earthquake was recorded by the seismic experiment in the Hangay Dome in scientific collaboration with the Lehigh University [15]. This network is comprised of 72 seismic broadband seismic stations located within and around the Hangay Dome from a latitude of ~44°N to ~50°N and a longitude of ~95°E to ~104°E ( Figure 2).

Figure 2. Distribution of the seismic stations network in Mongolia. Blue and red colored stars all permanent stations. Pink colored square is center of aimag. Green triangles are Hangay (Khangai) BB experiment temporary stations. Blue box represents study area along the South Hangay region
The stations recorded seismic data over a two-year period, June 2012 through April 2014. The dense spacing (~20-30 km between stations) geometry of array, as well as the length of data recordings, provide the opportunity to obtain better resolution of the subsurface.
The motivation of this study is to characterize seismicity and seismic source parameters in the South Hangay region. Highprecision earthquake locations and a complete, robust seismic catalogues are prerequisites for tectonic interpretation and seismic hazard assessment. Due to the intrinsically coupled hypocenter-velocity problem, earthquake locations and, in particular, their uncertainty assessment strongly depend on the velocity model. For reasons of faster and more robust ray tracing in 1D than in 3D velocity models [14] and also due to the lack of reliable 3D models because of inadequate coverage by high quality data in many regions, most seismological centers employ a priori 1D models for routine earthquake locations.

PMAS
Ideally, such 1D velocity models are obtained as solutions to the coupled hypocenter-velocity problem, by joint inversion for hypocenter and 1D velocity model by minimization of arrival time residuals [13].
In the Hangay region, different crustal thickness under the Hangay region have produced varying results depending on the methods used. Assuming the region is under isostatic equilibrium [16], the elevated tomography suggests a thickened crustal root. Seismic refraction measurements and Rayleigh wave phase velocities for the Hangay region indicate crustal thicknesses of 45-50 km with the lithospere possibly thinned to crustal thickness under Hangay [16]. A revised crustal thickness map by Zorin [17, Figure 8] based on consideration of seismic data, topography and gravity anomalies suggests that the South Hangay region is underlain by a crust up to 60 km thick with average values exceeding ~50 km [17,20]. In 2003, within the framework of the "Velocity structure of the Lithosphere on the 2003 Mongolian-Baikal Transect from SV waves" project, Mongolian scientists together with scientists from Russia and France determined crustal thickness as ~60 km [18].
Another study was carried out in the region in 2008. Results from the study, crustal thickness in region is determined as 48-50 km [19].
In recent studies, researchers are working on determining crustal thickness beneath each station using the receiver function method. From receiver functions, the average crustal thickness beneath recorded the seismic network in the study area is ~50 km. For this study, we assumed the crustal thickness is ~50 km for used 1D velocity model.
Consequently, we introduced 1-D model [11] with corresponding station corrections that may serve for routine uniform high-precision earthquake location [14] and as an initial reference model for 3D seismic tomography [12].

MATERIALS AND METHODS
The travel time of a seismic wave generated by an earthquake and recorded at a station is non-linearly dependent on hypocentral parameters and seismic velocities sampled along the ray path between source and station.
A minimum 1D velocity model with corresponding station corrections results from simultaneous inversions of a large number of travel times from selected high-quality events for both model and hypocenter parameters [11]. The method is designed to locate these events with the smallest possible uniform location error. The calculation of a minimum 1D model, following the routine procedure as defined in the Velest Users Guide [13] and Kissling et al., [12], is a trial and error procedure for different initial velocity models, initial hypocenter locations, and damping and control parameters for the coupled inverse problem. The term 'uniform location error' denotes that the sum of residuals for all events is minimized in the joint inversion. In a sense, the location accuracy is then relative to the full dataset [13]. For a set of well-locatable events with an azimuthal gap (the largest angular distance between two neighboring stations as seen from the epicenter) <180°, and at least 10 P-observations, which are evenly distributed within a network, the absolute location uncertainty can be estimated approximately by using randomly and systematically shifted hypocenters as initial locations and by analyzing the differences in final locations. In addition to these tests of hypocenter location accuracy, several tests have been conducted to assess the quality of the 1D velocity model. In my study region of South Hangay (45.7°-47.5°N, 97°-101°E), ~1700 events were selected by the experimental seismic stations network. Most of the recorded earthquakes are micro earthquakes while the largest earthquake occurred on 3 October, 2012 with local magnitude Ml=5.4 ( Figure 3). More than 500 aftershocks occurred in ensuing four months after this strong earthquake (Bayanbulag). The distribution of seismic stations and selected earthquakes for the data set are shown in Figure  3.

Figure 3. Seismicity map of South Hangay region for the period June 2012 -April 2014. Light blue triangular is Hangay BB seismic stations, pink star is permanent stations and red points are location of earthquakes which are scaled by magnitude
From the above events, 925 events were recorded in and around the South Hangay fault; 198 events with a magnitude greater than 1.5. We have selected events with Nsta>10 (Nstastation number), azimuthal GAP<160° (GAPthe largest angular distance between two neighboring stations as seen from the epicenter), RMS<0.4 (RMS-root mean square) from these earthquakes. Using these criteria, the dataset for inversion is composed of 120 earthquakes, with a total 9106 P and S-phase readings that uniformly distributed in the region of the South Hangay fault.
Calculation of a minimum 1D model for South Hangay Dome: We calculated a minimum 1D model, following the routine procedure as defined by Kissling et al., [12,13]. The data selection and the calculating process is described in Figure 4. In step 1, we determined an initial localization of 198 events by VELEST single-event mode scripts using the MNDC model. This MNDC is a half space model with a P wave velocity 6.11 km/s from surface to 45 km depth and 8.1 km/s below 45 km with a Vp/Vs of 1.73. The model was developed by Baljinnyam [21], which is still used today in routine seismic data analysis in MNDC. From the selected dataset for South Hangay fault system relocated with the updates priori 1D model (step 1, Figure 4).

Figure 4. Overview of procedure to calculate minimum 1D velocity model for South Hangay region (see text for details)
After that, we calculated the minimum 1D model using the selected dataset. Joint hypocenter determination calculation of several hundred events with VELEST is a valuable tool to identify errors in large travel time datasets [11]. The appropriate layering of the 1D model is estimated by a trial-and-error process. The priori model is calculated using an arbitrary layer thickness of 2 km for shallow crustal levels and increasing thickness with step 4 to 5 km to the Moho boundary which is fixed at 50 km depth ( Figure 5). Seismic wave travel time also depend on local geological conditions. Therefore, for the analysis we used station corrections for each station. In order to correct it qualitatively, we selected the station HD25 as a reference station, because it is located toward the center of the network and it provides a long high-quality record of data. This model and reference station are named the priori 1D model. We used 3 different starting velocity models with identical input and different control parameters to test the stability of the inversion process. A total 9106 arrival times of P and S phases from selected events were used in VELEST to determine a crustal velocity structure. With data from several stations available from a local or regional seismic source, the origin time can be determined by a very simple technique called Wadati diagram. Using this technique with criteria of RMS<0.4, Nsta>9, and R>0.8, we obtained an average Vp/Vs ratio of 1.73 with 0.03 standard deviation. The procedure started with 3 different initial models with average, higher and lower velocities with event locations, Proceedings of the Mongolian Academy of Sciences PMAS which were randomly distributed in their three spatial coordinates. The calculation of the minimum 1D velocity model comprised two inversion runs. In the first run, hypocentre locations were adjusted at every iteration, whereas seismic velocities and station delays were adjusted at every second iteration. Following Kissling [12], damping was set to 0.01 for hypocentre locations and station delays, and to 0.1 for seismic velocities. Figure  5a shows the set of initial velocity models for 1D inversion. Figure 5b represents the output velocity models from the various VELEST runs, after eight iterations. The inversion was stopped when model adjustments became insignificant (usually after 8 iterations). The goal of the first run was to find the appropriate minimum for each model in terms of seismic velocities, hypocentre locations, and station delays. After this iteration, minimum RMS model were selected for the next inversion run.

Figure 5. a) Various input initial velocity models for 1D inversion. b) Resulting velocity models after 8 iterations (starting initial model RMS=0.32, Low model RMS=0.36 and High model RMS=0.37). c) The output final model is marked the red line, updated a priori 1D model with corresponding station residuals shown as blue dashed line, black dashed line is our used MNDC velocity model
The second run, the calculation of the Minimum 1D model required multiple iterations to select and test control parameters appropriate to the dataset and to the problem. The damping factors provide a balance between the solution that minimize the errors and initial model. In this calculation, I considered an initial damping coefficient of 0.01 for hypocentral parameters, 0.1 for the station delays and 1.0 for the velocity parameters. Then, damping parameters for velocity variations and station corrections were selected optimizing the data of misfit function and the parameter resolution. The goal of this run is to minimize the total estimated location errors with fixed geometry.
We repeated the procedure for reduced number of layers by combining adjacent layers with approximate velocity values. In Figure 5c represents the P velocity models obtained at the end of the inversion procedure.
The range of possible velocity models obtained was tested with earthquake locations to select the best velocity model. This whole process leads to the final layering of a Minimum 1D model (Figure 5c, Figure 6 and Table 1). Obtained P wave 1D velocity model of the crust beneath the South Hangay region with a set of station correction is: A value of 6.16 km/s is obtained for the upper crust, and values around 6.27 and 6.63 km/s are obtained for the middle and lower crust, respectively ( Figure 6).

Figure 6. a) The obtained Minimum 1D P and S-wave velocity models for South Hangay region. b) Ray distribution in depth of 120 well-locatable events for inversion procedure. c) P and S-wave velocity ratio. The ratio between upper crust (1.76 down to 14 km depth) and lower crust (1.69)
The station corrections represent deviations of the velocity model. Positive and negative values correspond to local low and high velocity anomalies in the area of the recording station with respect to the station reference. In the present work, the station corrections were calculated for P and S waves. The values are in the interval of -0.7s to +0.5s for P waves and -0.8s to +0.6s for S waves (Figure 7).

Figure 7. Distribution of final station delays calculated with respect to the minimum 1D model for South Hangay region. Colored open circles indicate the positive and negative values of delays, respectively, for stations with at least 10 P-wave observations, relative to the HD25 reference station that is marked by a pink star
In global view, the northern part of the network is dominated by positive corrections, indicating true velocities lower than the model velocities. The South Hangay region in general has negative corrections, here the true velocities are faster than the model velocities.
In this work, we did not consider correlation between final station correction and station elevation, because in the 1D inversion we considered station altitude.

RESULTS AND DISCUSSION
We performed a series of simultaneous inversion using VELEST and the arrival time data obtained by MNDC-IAG, in order to determine P and S wave 1-D velocity structures of the crust and upper mantel beneath the South Hangay region and with a set of station corrections. The computed 1D model is shown in Figure 8 (Output 1D inversion) and we compared the model with some international standard models: IAG (Institute of Astronomy and Geophysics), IASPEI91 (International Association of Seismology and Physics of the Earth's Interior), RSTT (Regional Seismic Travel Time), IIEES (International Institute of Earthquake Engineering and Seismology).

P-wave station correction S-wave station correction
Proceedings of the Mongolian Academy of Sciences

Figure 8. The obtained Minimum 1D velocity model compared with International standard models
The obtained 1D velocity model shows similar velocity gradients with IIEES velocity model at 20-30 km depth. From 30-35 km depth, the model shows similar velocity gradients with IASPEI91 model (Figure 8). From the model, we can see an intriguing feature at at a depth of 45 km. P velocity value shows 7.16 km/s at this depth. This feature may be explained by a third, basalt crustal layer. In many stable continental interiors there is a third, basalt crustal layer with a velocity of 6.8-7.2 km/s. The seismic velocity below Moho (Pn velocity) is typically about 8 km/s [22].
Using this inversion result, I relocated hypocenter of 1700 selected events by program HYPOCENTER provided by the SEISAN software package. The obtained results were compared with the hypocenter location using the current velocity model which is being used for the routine data analysis by the Mongolian network.
As mentioned before, it shows clearly that the depth is not constrained, 80% of event's depth were fixed at 2 km in the NDC catalogue. After the hypocentre localization procedure, we significantly improved the hypocentres with a considerable reduction of depth uncertainty (Figure 9). The distribution of seismicity confirms that the majority of the seismic activity is clustered in space in the South Hangay fault area. The estimated 1-D velocity model showed a reduction of total RMS value of all events calculation. Relocated and initial hypocenters are plotted in Figure 9. The epicentral coordinates of the events are well resolved with error less than 4.0 km on average. The largest shifts between initial and relocated hypocenters were observed at the south west and north-west side of the map, which is a total of 4 percent from all events. These large shifted locations might related to sparse seismic network and small events that difficult to fix seismic phases ( Figure 10).

Figure 10. Comparison between initial and relocated hypocenters. Blue plus symbols mark initial hypocenters, red symbols are relocated hypocenters. Black lines are active tectonic faults, pink line is border of provinces, blue triangles are Hangay seismic experiment stations. The histogram shows original and relocated epicenter difference
The new 1D minimum model presented in this paper can be used in a straightforward way to locate earthquake in region. Considering our results, the use of a greater number of layers seems the be unjustified. Finally, purposed this model for future tomographic studies and seismic hazard assessment in this region.

CONCLUSIONS
We performed a series of simultaneous inversion to estimate a 1D velocity model with station corrections of the South Hangay region in central Mongolia. For this purpose, we used an efficient method for estimating the minimum 1-D velocity model by solving the coupled hypocenter-velocity model problem.
The method is minimized P-wave and S-wave residuals of the data set according to the procedures outlined by Kissling et al.,[9] that inverted a set of earthquake data by VELEST program.
The simultaneous inversion for structure and hypocenters was carried out using 198 detected events with local magnitude more than 1.5 from over 925 events which occurred along in and around the South Hangay fault system that recorded by the BB experiment.
Consequently, we obtained P wave velocity 6.16 km/s, 6.27 and 6.63 km/s for the upper, middle and lower crust, respectively.
The results are compared with several existing local and global model values.
Using our proposed model, we relocated the ~1700 selected events and compared MNDC locations. The result shows that the use of obtained model is improving location and depth of the events.
Therefore, we proposed this new model for routine hypocenter determination for this region, which is expected to improve the accuracy of the locations and act as an initial reference model for seismic tomography study.