Dual threshold search method for asperity boundary determination based on geodetic and seismic catalog data

Xiohng Wng , Zhongzheng Zhou , Cijun Xu , Yngmo Wen , Hu Liu

a School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China

b Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China

c Key Laboratory of Geophysical Geodesy, Ministry of Natural Resources, Wuhan University, Wuhan 430079, China

Keywords:GPS Earthquake catalog Dual threshold search method Asperities Haiyuan fault

ABSTRACT As an important model for explaining the seismic rupture mode, the asperity model plays an important role in studying the stress accumulation of faults and the location of earthquake initiation.Taking Qilian-Haiyuan fault as an example, this paper combines geodetic method and b-value method to propose a multi-source observation data fusion detection method that accurately determines the asperity boundary named dual threshold search method. The method is based on the criterion that the b-value asperity boundary should be most consistent with the slip deficit rate asperity boundary. Then the optimal threshold combination of slip deficit rate and b-value is obtained through threshold search,which can be used to determine the boundary of the asperity.Based on this method,the study finds that there are four potential asperities on the Qilian-Haiyuan fault: two asperities (A1 and A2) are on the Tuolaishan segment and the other two asperities (B and C)are on Lenglongling segment and Jinqianghe segment, respectively. Among them, the lengths of asperities A1 and A2 on Tuolaishan segment are 17.0 km and 64.8 km, respectively. And the lower boundaries are 5.5 km and 15.5 km, respectively; The length of asperity B on Lenglongling segment is 70.7 km,and the lower boundary is 10.2 km.The length of asperity C on Jinqianghe segment is 42.3 km, and the lower boundary is 8.3 km.

1. Introduction

In order to describe the rupture mode of large earthquakes and main shocks on active faults, Kanamori and Stewart [1] first introduced the concept of asperities into the field of seismology when studying the Guatemala earthquake in 1976. They believed that asperities are areas with high stress on active faults.And these areas have the properties of high roughness, high friction and low pore pressure. For seismic risk assessment, the location of the asperities is usually the starting point of potential earthquake rupture[2].For the asperities that have experienced large earthquakes,they can be determined by seismic inversion and dislocation model.However, on the active fault segmentations where no major earthquakes have been observed, how to identify the asperities is still a hot issue of people's attention.

There are two main methods to determine the distribution of asperities on faults during the inter-seismic period.The first is the geodetic based slip deficit method (SD method), according to the definition of asperity,it is generally considered that the area where the slip deficit is concentrated on the fault is the asperity.Geodetic data can be employed to obtain the concentrated area of slip deficit on the fault to determine the distribution of the asperity[3-5].The second method is to use b-value to identify asperities(BV method).Based on the inverse relationship between b-value and stress[6,7],Wiemer and Wyss[8].proposed to use the abnormally low b-value calculated from the background seismic activity to outline and determine the location of the asperities.This method has been used to analyze the risk of strong earthquakes in different segments of active fault zones [9-11].

Different methods for determining asperities have significant differences in the discrimination method, data processing and criteria. Combining different methods can effectively improve the reliability of estimated asperities. This research proposes a dual threshold search method to determine the distribution of asperities by combining the SD method and the BV method. By minimizing the difference of the asperity boundaries, the asperity boundary determined by BV method is used to constrain the optimal asperity distribution threshold based on the result of SD method.Compared with single SD method or BV method, this method not only effectively improves the reliability of the determination of asperities,but also gives an optimal threshold for determining the range of asperities. This research will take Qilian-Haiyuan fault as an example;the accurate distribution of asperities of the fault will be given based on the dual threshold search method.

2. Study area and data

As one of the most intense seismic activity and tectonic deformation zone on the northeastern margin of the Qinghai-Tibet Plateau, the Qilian-Haiyuan fault is surrounded by the Alxa block,the Ordos block, the Lanzhou block and the Qilian block (Fig. 1).Previous studies have shown that the Qilian-Haiyuan fault is a large Left-lateral strike-slip fault [12], and the east and west sides are respectively connected with the Liupanshan fault and the western margin of Qilianshan fault. The Haiyuan fault includes the westmiddle-east segments. The Qilian-Haiyuan fault includes the Tuolaishan fault, Lenglongling fault, Jinqianghe fault, Maomaoshan fault and Laohushan fault. This research took the Qilian-Haiyuan fault as the target fault. In 1920, the M8.5 Haiyuan earthquake occurred on the eastern part of the Qilian-Haiyuan fault zone,which caused a surface rupture of 240 km and a maximum left-lateral coseismic displacement of 10 m [13]. After that, the 1927 M8.0 Gulang earthquake occurred near the Qilian-Haiyuan fault zone[14].Since then,the Qilian-Haiyuan fault zone has been in a stable state,and there is a 220 km long potential seismic hazard zone: the Tianzhu earthquake gap [13] between the two major earthquakes.Until January 8, 2022, a strong earthquake (M6.9) shocked the-Menyuan County in Qinghai Province (Fig.1). This earthquake was located near the junction of the Tuolaishan segment and the Lenglongling segment(U.S.Geological Survey,USGS).

The data used in this study includes GPS data and earthquake catalog data. Among them, the GPS data comes from Wang et al.[15] and we selected the GPS data in the range of longitude 93.2E -113.8E and latitude 32.2N -44.6N, and finally got 586 stations(Fig.2).The seismic catalog data contains seismic data from the 1920 Haiyuan earthquake to May 2021.Fig.1 shows the seismic events with M ＞2.0.

3. Dual threshold search method

Slip deficitrepresents the rate of stress accumulation on the fault plane,the greater rate of deficit,the faster accumulation of stress on the fault plane.According to the definition of asperity,it is generally considered that the area where the slip deficit concentrates on the fault zone is the asperity. The main models for obtaining the interseismic slip deficit from geodetic data are:2D model and 3D model.Compared with the 2D inter-seismic model, the 3D model has the advantages of variable fault strike,consideration of the“translation”and“rotation”of the block,and the influence of secondary faults in the modeling process, and it has been widely used in many fault zones [5,16-18].

Fig. 1. The regional structure of the Qilian-Haiyuan fault and the distribution of historical earthquakes. TLS: Tuolaishan segment; LLL: Lenglongling segment; JQH: Jinqianghe segment; MMS: Maomaoshansegment; LHS: Laohushan segment; HYF: Haiyuan segment; LPS: Liupanshan segment. The blue line is the Tianzhu earthquake gap.

Fig.2. GPS velocity field on the northeastern margin of the Tibet Plateau.(a)The total GPS velocity field of the study area,relative to the Eurasian plate.A:Alxa block,B:Qilian block,C: Lanzhou block, D: Ordos block; (b) The velocity field around the Qilian-Haiyuan fault, relative to the Eurasian plate. The red line is the location of the Qilian-Haiyuan fault.

The b value is a measure of seismic activity,and it can represent the stress state of the study area. Lower b value means stronger seismic activity,and higher stress state.According to the hypothesis of the asperity model, the generation and occurrence of strong earthquakes are generally located on locked areas or asperities with high stress accumulation [8,19]. Therefore, the changes in the spatial distribution of b-value have become a means to effectively assess the mid-to-long-term seismic risk [20,21], which provides guidance for the division of the spatial distribution of asperities.

Fig. 3 is a flowchart of the dual threshold search method. This method is mainly divided into 5 steps:First,obtain the distribution of slip deficit on the fault plane and the distribution of b value along the surface fault line according to the geodetic data and the seismic catalog data,respectively.Secondly,set the initial inter-seismic slip deficit rate threshold and b-value threshold to ensure that the initial threshold can completely cover the distribution range of asperities. For the b-value, the initial thresholdsare based on the local maximum. For the slip deficit, simulation experiments are required to determine the initial threshold. In this study, we took 30%of the maximum slip deficit as the initial threshold which can completely cover the range of asperities. Next, based on the initial slip deficit threshold and the initial b-value threshold, expanding the search ranges [Sd, Sd] and [bm, bm], and give the corresponding search step size ΔSd and Δbm, respectively. Then,through iterative calculation of different combinations of Sd and bm in the search range, determine the optimal threshold combination (Sd, bm) based on the most consistent boundary of the asperities from slip deficit and b-value. In the process, the criterion is realized by minimizing the sum of squared differences between the east-west boundary (Eand W, Eand W) of the asperities from SD method and BV method. Finally, the optimal threshold combination (Sd, bm) is used to determine the final boundary of the asperities.

4. Determination of the asperity distribution on Qilian-Haiyuan fault

4.1. Asperity distribution on Qilian-Haiyuan fault based on slip deficit model

Based on the GPS data, this study employed a linear spherical block model [3] to calculate the distribution of slip deficit on the Qilian-Haiyuan fault. The linear block theory, which decomposes surface velocity fields into four components:(1) plate rotations, (2)elastic deformation from faults with kinematically consistent slip rates, (3) elastic deformation from faults with spatially variable coupling, and (4) homogeneous intrablock strain. Elastic deformation rates are computed for each fault segment in a homogeneous elastic half-space using multiple optimal planar Cartesian coordinate systems to minimize areal distortion and triangular dislocation elements to accurately represent complex fault system geometry.Block motions, fault-slip rates, elastic coupling, and internal block strain rates are determined simultaneously using a linear estimator with constraints from both geodetically determined velocity fields and geologic fault-slip rate estimates. According to previous studies[22,23], this study divided the northeastern margin of the Tibet Plateau into Ordos block,Alxa block,Qilian block and Lanzhou block(Fig.2(a)).Since most of the GPS stations are far away from the fault,and the GPS data is not sensitive to the difference in fault geometry at depth.Therefore,we set the dip angle of the Qilian-Haiyuan fault to 65[24], and assume other faults in the model to be vertical[15,25-28]. At the same time, the locking depth was restricted by fitting GPS data [29,30]. In this study, the average locking depth of the major faults on the northeastern margin of the Tibet Plateau was first determined to be 20 km(Fig.4(a)),and then the optimal locking depth of the major faults in the study area was determined based on this depth (Fig. 4(b)-(h)), which is consistent with the result given by Li et al.[29].

Fig. 3. Flow chart of dual threshold search method. The method consists of three parts, namely Data and Model, Parameter Calculation and Threshold Determination.

Then, we inverted the slip deficit on Qilian-Haiyuan fault. To ensure that the depth of the model is below the locking depth of the Qilian-Haiyuan fault, the depth of the fault model was set to 24 km[22,24].In this study,the optimal smoothing factor(β)was determined based on the slip rate of the Tuolaishan fault, the maximum, minimum and average slip rates of this fault were calculated by varying the value of β. Then, the β that makes the slip rates of this fault comparable to the rate obtained from previous block model inversion was considered reasonable [28].As shown in Fig. 5, when β = 1, the maximum, minimum and average changes of the slip rate of Tuolaishan fault tend to be flat and close to the slip rate(dashed green line in Fig.5(a))estimated from the previous block motion, so β=1 was determined as the optimal smoothing factor.

Fig. 4. Optimal locking depths of main faults on the northeastern margin of Tibet plateau. (a) The average locking depth on the northeastern margin of the Tibet Plateau; (b) the locking depth of the Tuolaishan fault;(c)the locking depth of the Lenglongling fault;(d)the locking depth of the Jinqianghe fault;(e)the locking depth of the Maomaoshan fault;(F)the locking depth of the Laohushan fault; (g) the locking depth of the Haiyuan fault; (h) the locking depth of the Liupanshan Fault.

Geodetic inversion (Fig. 6) shows that there are mainly five obvious slip deficit areas on Qilian-Haiyuan fault, two asperities are on the Tuolaishan segment and the other three asperities are on Lenglongling segment, Jinqianghe segment and Laohushan segment,respectively.The maximum slip deficit rate is 7.8 mm/yr in Tuolaishan segment. Based on the dual threshold search method, we first delineated the boundary of the asperity with a smaller slip deficit threshold (Fig. 6), with the model and smoothing factor determined, we had done some simulation experiments and found that for a single asperity, taking 30% of the maximum slip deficit as the initial threshold can completely cover the range of asperities.The initial slip deficit threshold Sd was set to 2.3 mm/yr and 2.4 mm/yr for the two asperities on Tuolaishan segment and 2.3 mm/yr, 2.1 mm/yr and 2.0 mm/yr for the other three asperities on Lenglongling segment,Jinqianghe segment and Laohushan segment, respectively. It is worth noting that these thresholds are only preliminary result, the precise threshold for delineating the asperity boundary will be obtained from the dual threshold search method.

Fig. 5. Results of the smoothing iteration test. (a) A full block model was run using a variety of values of β,with all other parameters held constant.The extrema are shown as black lines, and the mean is shown as a red line. The relative blockmotion velocity calculated on the Tuolaishan fault is shown by thehorizontal dashed green line.(b)Plot of statistics of the modeled velocity field as a function of β.The blue line represents the mean residual velocity magnitude (speed), and the orange line shows the percentage of residual speeds less than the meanobserved uncertainty speed.

Fig.7(a)shows the observation and simulation of the GPS data,Fig. 7(b) represents the residuals which indicates that the simulation and observation are in good agreement, and about 82% of the GPS data have residuals less than 1 mm/yr.

4.2. Asperity distribution on Qilian-Haiyuan fault based on b-value

In this study, the linear regression method was employed to calculate the b-values, considering the large dip of the Qilian-Haiyuan Fault, the b-values in planar space have limited constraints of the asperities on the fault plane. Therefore, only the bvalues on the surface fault line were calculated to constrain the asperities, we took the sampling points at an interval of 20 km on the fault line to calculate the b value and the size of space grid is the area of the circle with a radius of 10 km centered on the sampling points on the fault line.Two large earthquakes have been recorded near the Qilian-Haiyuan fault zone(Fig.1),the most recent of which was the M8.5 Haiyuan earthquake in 1920 [31] and the M8.0 Gulang earthquake in 1927 [14]. The occurrence of large earthquakes releases stress accumulation during the seismic cycle,which has a large impact on the distribution of seismicity in the study area and also affects the results of b-values. Since the seismogenic faults of the 1927 Gulang earthquake were not on the Qilian-Haiyuan fault, and the Gulang event was far from these segments, the b-values on the Tuolaishan segment to the western part of the Maomaoshan segment were calculated based on the catalog from the record to 2021.However,there is a serious lack of earthquake catalogues on the eastern part of Maomaoshan segment to Haiyuan segment. In most areas, the number of earthquakes within 10 km from the sampling point is less than 30,the b-value results are less reliable.Therefore,b-value calculation is not performed for this part, and the asperities on this part is determined by geodetic inversion. Due to the relatively small amount of seismic data in this area, we used the M(minimum magnitude of completeness) [32] of the fault segments to replace the Mof the sampling points on the segment in this study,and the Mwas set to 2.0, 2.1,1.8 and 1.9 on Tuolaishan, Lenglongling, Jinqianghe and Maomaoshan segments, respectively(Fig. 8).

Fig. 6. Slip deficit distribution and classification of asperities in Qilian-Haiyuan fault based on slip deficit.

The result indicates that there are two low b-value regions (A1 and A2)on the east and west sides of the Tuolaishan segment(Fig.9(a)),here we set the initial threshold of the two asperities to 0.755(A1)and 0.819(A2)based on the local maximum b-value.There is also an obvious low b-value region (B) on Lenglongling segment and the initial threshold was set to 0.816.The b-value drops rapidly on the middle of the Jinqianghe segment and reaches its lowest point on the middle part of the segment(Fig. 9(b)). Therefore, we can infer that there may be an asperity (C) on Jinqianghe segment and we set the initial threshold to 1.032 based on the local maximum b-value. Since there are no b-value results on the Maomaoshan segment to Haiyuan segment, we did not discuss the results of this part in this study. The thresholds for the division of asperities in each segment of the Qilian-Haiyuan Fault based on bvalue are shown in Table 1.

4.3. Determination of asperity distribution on Qilian-Haiyuan fault via dual threshold search method

Fig.7. Modeled and residual GPS velocities of the northeastern margin of Tibet Plateau.(a)GPS observations and simulations;(b)GPS residuals,the attached picture is the residual statistics. A:TLS, B:LLL, C:JQH, D:MMS, E:LHS, F:HYF.

Fig. 8. Magnitude distributions and GR b-value estimates. Minimum magnitudes are marked by the vertical dashed lines. (a) Tuolaishan segment; (b) Lenglongling segment; (c)Jinqianghe segment; (d) Maomaoshan segment.

Fig.9. Distribution of asperities in each segment of the Qilian-Haiyuan fault based on b-value. (a) Tuolaishan-Lenglongling segment;(b)Jinqianghe-Laohushan segment;(c)Haiyuan segment.

Table 1 B-value-based classification threshold of asperities in each segment of Qilian-Haiyuan fault.

Aftertheinitialrangeofasperitieswasdetermined,thisstudy used the dual threshold search method to accurately divide the asperities on each segment (Fig. 10). By comparing the results between slip deficit model and b-value, it can be found that asperities can be determined by both two methods on Tuolaishan, Lenglongling and Jinqianghe segments. However, the BV method cannot determine whether there are asperities on the Laohushan segment due to the lack of seismic catalogues. Therefore, this study only employed the dual threshold method to determine the distribution of the asperities onTuolaishan,LenglonglingandJinqianghesegments.Asaresult,two asperities(A1 and A2)were determinted on the Tuolaishan segment and the other two asperities (B and C) were determined on Lenglongling segment and Jinqianghe segment,respectively.Among them,the lengths of asperities A1 and A2 on Tuolaishan segment are 17.0 km and 64.8 km, respectively. And the lower boundaries are 5.5 km and 15.5 km, respectively; The length of asperity B on Lenglongling segment is 70.7 km,and the lower boundary is 10.2 km.The length of asperity C on Jinqianghe segment is 42.3 km,and the lower boundary is 8.3 km. The thresholds for the division of asperities in each segment of the Qilian-Haiyuan Fault based on dual threshold search method are shown in Table 2.

5. Discussio n

To ensure the accuracy of the slip deficit model inversion result,a checkerboard test was played in this study. Based on the inversion result,we set up checkerboards of 40 km×12 km,40 km×10 km,and 40 km×5 km in this experiment.At the same time,different slip deficit rates were set for analysis based on the asperities on the Qilian-Haiyuan fault,and the smoothing factor(β)was set to 1.It was found that the difference in slip deficit rate between 6 mm/yr and 8 mm/yr had no effect on the checkerboard test results. On the premise that the length of checkerboard is 40 km and the width is greater than 12 km, the GPS data can well constrain the slip deficit model(Fig.11(a-b)).When the width of checkerboard is 10 km,the GPS data can only well constrain the slip deficit within 10 km of the fault in depth(Fig.11(c-d)).When the checkerboard decreases again,the GPS data can only well constrain some areasin the shallow region,such as the western part of Tuolaishan segment and Jinqianghe-Laohushan segment (Fig.11(e-f)). In brief, the GPS data employed in this study can well constrain the slip deficit model(Fig.6).

Previous studies did not give the precise boundaries of the asperities on the Qilian-Haiyuan fault [24,33-35]. In this study, the optimal threshold combination of b-value and interseismic slip deficit rate was obtained based on the dual threshold search method, and the specific location and boundaries of the asperities on the Qilian-Haiyuan Fault were given.The results show that there were four asperities on the Qilian-Haiyuan Fault. The slip deficit areas are mainly concentrated on the west side of the Qilian-Haiyuan fault. And it is worth noting that the slip deficit near the epicenter of the 2022 Menyuan earthquake is obvious, which is consistent with the results of previous studies [24,33-35]. It is worth noting that Cavaliˊe et al. [33] and Jolivet et al. [34,35]employed InSAR data to study the Qilian-Haiyuan fault and found a 35 km long creep zone in the“Tianzhu Earthquake Gap”,but we did not detect a shallow creep in this area. Considering that the resolution of InSAR observations is much higher than GPS observations,and as shown in Fig.2,GPS stations are sparsely distributed in the near field of the Qilian-Haiyuan Fault.Therefore,the reason for this difference may be due to the insufficient data resolutionof the GPS observations.

Fig.10. The distribution result of asperities on Qilian-Haiyuan fault based on dual threshold search method.

Table 2 The details of asperities in each segment of Qilian-Haiyuan fault based on dual threshold search method.

Fig.11. Checkboard test with checkerboards of 40 km × 12 km, 40 km × 10 km, and 40 km × 5 km and different slip deficit rates.

Different from other segments of the Qilian-Haiyuan fault,the geodetic inversion indicates that there is no significant slip deficit on the Haiyuan fault, only weak slip deficit exists in the shallow depth(about 5 km).There is no slip deficit in the deep part,which is close to a creeping state (Fig. 6), indicating that the stress accumulation of this region is slower than that of other segments.Since the Haiyuan fault is the seismogenic fault of the 1920 Haiyuan earthquake[13,36,37],the slower stress accumulation in this segment may be due to the influence of the Haiyuan earthquake.The stress accumulated in the fault was released during the Haiyuan event.Therefore,the seismic risk of this segment is low at present.

6. Conclusion

This study combines the geodetic method and b-value method to propose a multi-source observation data fusion detection method that accurately determines the asperity boundary named dual threshold search method.The method is based on the criterion that the b-value based asperity boundary should be most consistent with the slip deficit rate asperity boundary. Then the optimal threshold combination of slip deficit rate and b-value is obtained through threshold search, which can be used to determine the boundary of the asperity. Taking the Qilian-Haiyuan fault as an example,the dual threshold search method was used to determine the precise boundary of the asperities,and it is found that there are four potential asperities on the Qilian-Haiyuan fault:two asperities(A1 and A2) are on the Tuolaishan segment and the other two asperities (B and C) are on Lenglongling segment and Jinqianghe segment, respectively. Among them, the lengths of asperities A1 and A2 on Tuolaishan segment are 17.0 km and 64.8 km, respectively. And the lower boundaries are 5.5 km and 15.5 km, respectively;The length of asperity B on Lenglongling segment is 70.7 km,and the lower boundary is 10.2 km. The length of asperity C on Jinqianghe segment is 42.3 km,and the lower boundary is 8.3 km.There are no potential asperities on the Haiyuan segment, which may be due to the fact that the 1920 Haiyuan earthquake has released most of the accumulated stress.

Author contributions

Xiaohang Wang, Zhongzheng Zhou, Hu Liu performed the experiments, Xiaohang Wang wrote and revised the manuscript.Caijun Xu and Yangmao Wen designed the study, analyzed the experimental results and revised the manuscript.

Con

fl

icts of interest

The authors declare that there is no conflicts of interest.

Acknowledgments

This work is supported by the National Key Research and Development Plan of China under Grants No.2018YFC1503604,the National Natural Science Foundation of China under Grants No.41721003, No.42074007, the Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University,No.19-01-08. Some figures in this paper were prepared using Generic Mapping Tools [38].