Quantitative Study of Crustal Structure Spatial Variation Based on Gravity Anomalies in the Eastern Tibetan Plateau: Implication for Earthquake Susceptibility Assessment

Areas with sharp spatial changes in crustal structure are usually considered prone to strong earthquakes, especially late Cenozoic structural zones. However, quantitative relationships have not been established between the occurrence of earthquakes and the extent to which the crustal structure changes. The crustal structure can be reflected by gravity anomalies. In this study, we investigate the crustal structure variation in the eastern Tibetan Plateau based on the gravity data acquired by the Gravity field and steady‐state Ocean Circulation Explorer satellite. The multiscales wavelet decomposition and power spectrum methods are used to analyze the spatial variation of crustal structure at different depths. Based on a fourth‐order wavelet map, four variable factors are constructed to represent the spatial variations of the crustal structure, and their statistical relationship with earthquakes is studied. The results show that the changing rate of the gravity anomalies is significantly positively correlated with the occurrence of earthquakes. Areas of high earthquake hazards are places where the gravity field changes dramatically, usually with a value of 11°–18° for the SLP (slope) index. The assessment efficiency of the SLP increases with the earthquake magnitude. According to the SLP index, several regions are assessed to have a high level of risk for strong earthquakes in the eastern Tibetan Plateau, including the Longmenshan Fault belt, where a Ms8.0 earthquake occurred in 2008, thus demonstrating the promising implications of our method in the earthquake susceptibility assessment.


Introduction
Geophysical data play an important role in seismic risk assessments and predictions for providing information about the geological setting of the deep crust. Gravity is one of the most important types of data among the kinds of geophysical data. The heterogeneity of crustal density leads to spatial variation of the gravity field. Bouguer gravity anomalies can reflect density changes in the crust after removing the topographic mass effect and applying the free-air reduction (Heiskanen & Moritz, 1967). The density variations in the crust are often closely related to faults and tectonic structures, which provide geological setting information about seismic risk in research.
By applying wavelet multiscale decomposition (Mallat & Hwang, 1992;Mallat & Zhong, 1992) and the power spectrum method (Spector & Grant, 1970;Zhang & Zhao, 2007), the Bouguer gravity anomaly is usually used to study the active fault distributions at different depths and the seismicity of active earthquake areas Jiang et al., 2014;Tian et al., 2017) and to calculate the Moho depth (Aitken et al., 2013;Meliani et al., 2016;Moritz, 1990;Shin et al., 2015;Sjöberg, 2009Sjöberg, , 2013Sjöberg & Bagherbandi, 2011;Sobh et  This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. et al., 2012). Earthquakes occur as a result of intense fault activity; therefore, earthquakes often occur along active faults. Fault length and penetration depth are important factors affecting earthquake magnitude. The crustal structure changes, especially the crustal thickness change, indicate the scale size of the faults along the changing belt. For this reason, the gradient zone of the Moho depth is usually considered prone to strong earthquakes Li, Fang, & Braitenberg, 2017;Liu et al., 2016;Wei et al., 2016;Wu et al., 2015), which is also true for the zone of the crustal density transition area indicated by gravity anomalies (Yang et al., 2012;Li, Fang, & Braitenberg, 2017) and the Curie depth gradient belt (Gao et al., 2015). A numerical simulation of the generation mechanism of large earthquakes shows that the Moho upheaval is coincident with the high stress concentration rate locations (Liu et al., 2016).
However, the relationship between the spatial distribution of strong earthquakes and the crustal structure spatial variation is mostly qualitatively described based on map observations. The quantitative description of these relationships can provide more insight for research on seismogenic geology settings. In this paper, we study the crustal structure characteristics of the eastern Tibetan Plateau at different depths based on the gravity data calculated from the GECO model (Gilardoni et al., 2016), which combine the signals from EGM2008 and the GOCE satellite data in low-medium degrees. Furthermore, a quantitative analysis method of the seismogenic geology setting study is proposed based on the gravity data that can also be used for numerical research of the seismic tectonic backgrounds of other geophysical fields, such as the crustal velocity structure.

Geological Setting
The Qinghai-Tibet Plateau (QTP) is a young plateau that has undergone large-scale and high uplift during the collision and convergence of the Indian and Eurasian plates since 65-50 Ma (Molnar & Tapponnier, 1975;Harrison et al., 1992;Molnar et al., 1993;Tapponnier et al., 2001;). The QTP and its surrounding areas show strong three-dimensional changes in terms of the crustal thickness, seismic wave velocity, resistivity, and topography (Bai et al., 2010;Chen & Molnar, 1981;Lu et al., 2019). The spatial changes in the crustal structure are also reflected in gravity anomaly data. Because of the uplift of the QTP, the crust under the main part of the plateau has thickened, and the bottom of the crust is sunken, which produces an area with very significant negative Bouguer gravity anomalies ( Figure 3) and high-elevation topography (Figure 1). Because of the strong tectonic activity, the crustal structure shows apparent lateral differentiation in this area and is divided into multiple tectonic units by several active faults ( Figure 1). The seismic frequency is very high in the QTP and its periphery due to the variety of tectonic types and tectonic activity in this area.
Along the eastern boundary of the QTP, the Longmenshan fault belt (LFB) divides the QTP and the Sichuan Basin (SB); the Wenchuan Ms8.0 earthquake occurred on this belt in 2008. Due to the uplift of the QTP and the eastward crustal material movements escape blocked by the SB, the relief across the LFB changes dramatically. West of the LFB, the upper crust of the Longmenshan subblock (LSB) was severely destroyed and the lower crust was thickened during the uplift of the QTP , and it shows obvious magnetic, Curie depth, and velocity differentiations from the crust of the SB (Gao et al., 2015;. The Bayan Hara block is in the western part of the LSB, and the Longriba fault between them is a distinct velocity gradient belt (Zhang & Wei, 2011;Lei et al., 2019;). Broadband magnetotelluric data show that the resistivities of the SB-LSB and the West Qinling orogenic belt are very different (Xue et al., 2019).
The Sichuan-Yunnan block (SYB) is a rhomboidal block with strong seismic activity that was extruded during the uplift of the QTP and is confined by the Xianshuihe fault, the Anninghe fault, the Zemuhe fault, the Xiaojiang fault, and the Honghe fault. The Lijiang-Xiaojinhe fault (LXF) divides the SYB into two subblocks named the Northwest Sichuan subblock (NSSB) and the middle Yunnan subblock (MYSB). The maximum principal stress directions in the SYB show a clockwise rotation from east-west to southeast and southwest (Li, Liu, et al., 2017;Zhao et al., 2012). After the Ms8.0 Wenchuan Earthquake, the slip along the LXF changed from a left-lateral slip to thrust shortening (Zhao et al., 2012). The low-velocity crustal flow area is confined by the LXF . The crustal structure shows significant differences between the two subblocks in the SYB.
The Ordos block (OB) is surrounded by a series of narrow grabens Liang et al., 2018;Zhang et al., 2017) and shows a nearly homogeneous crustal structure inside the block, with an average crustal thickness of 41.5 km. In the north margin, the Yinchuan-Hetal graben is located in the Paleoproterozoic khondalite zone, with a high Poisson's ratio and a thickened crust. In the west margin, the crustal thickness and Poisson's ratio are changed significantly across the Liupanshan thrust belt Wang et al., 2017).

Data
The earthquake catalog used in this paper is provided by the China Earthquake Data Center (http://data. earthquake.cn). It covers from an area from 88°E to 110°E and 20°N to 44°N, a time period from 1970.01 to 2019.06, and magnitudes greater than 3.0. The earthquakes outside the intersection coverage of the active faults and gravity data are deleted to implement further comprehensive statistical analyses. The aftershocks are removed with the G-K method (Gardner & Knopoff, 1974).
The gravity data grids used in this paper are the Bouguer gravity data downloaded from the International Center for Global Earth Models (http://icgem.gfz-potsdam.de/calcgrid) with the GECO model (Gilardoni et al., 2016) and gravity_anomaly_bg selected in tide-free mode. The (simple) Bouguer gravity anomaly is defined as the classic gravity anomaly minus the attraction of the Bouguer plate. The topographic heights H(λ,φ) are calculated from the ETOPO1 model. These data cover 89°E to 109°E and 21°N to 43°N spatially with a resolution of 10 km. After the wavelet decomposition analysis, the raster data outside the intersections area of the active faults were extracted as further statistical data.
The main active fault lines used in this study are simplified from the latest active fault data provided by the Active Fault Survey Data Centre at the Institute of Geology, China Earthquake Administration (Xu et al., 2016). Combining the coverages of active fault data and gravity data, the final specific research area of this study is shown in Figure 3.

Multiple Scale Wavelet Decomposition Method
The wavelet decomposition method (Mallat & Hwang, 1992;Mallat & Zhong, 1992) was used to obtain different scales of the gravity field source signals.
For the two-dimensional gravity anomaly field, where The term A j f(x,y) represents the low-frequency approximation part of the two-dimensional multiscale analysis space, and the terms D h j f x; y ð Þ, D v j f x; y ð Þ and D d j f x; y ð Þ represent the high-frequency details in the horizontal, vertical, and diagonal directions, respectively. The function ∅ is the scaling function, and ψ is the wavelet function. The terms m 1 and m 2 represent the discrete coordinate values of two dimensions in the x and y directions, respectively, and Z represents a set of integers. Using the gravity anomaly data on the grids as the coefficients {s 0 } of the spatial scale function, we can calculate the coefficients {s j } of the other spatial scale functions and the coefficients {d h,j }, {d v,j }, and {d d,j } of the wavelet function to realize the multiscale decomposition of 2-D gravity anomalies. This decomposition is simply written as follows: where A 4 G is the fourth-order approximation part of the wavelet decomposition. D j G is the sum of the three detail parts, namely, horizontal, vertical, and diagonal, and the high-frequency detail part of the jth-order wavelet decomposition. The anomalies on the low-order wavelet detail map represent the density anomalies in the shallow layer of the crustal, while the high-order details are opposite. The average depth of the field source can be estimated by the power spectrum method.

Power Spectrum Method and the Average Field Source Estimate
The average depth of field source of detail maps derived from the wavelet decomposition method can be estimated using the power spectrum method (Spector et al., 1970). The average depth can be written as follows: where Z(ρ) is the power spectrum and ρ is the wave number within 1 km. The power spectra and the average depths of the firstto fourth-order wavelet detail map (1-D-4-D) are shown in Figure 2.

The Bouguer Gravity Anomaly
The eastern QTP shows low-negative Bouguer gravity anomalies, while the tectonic units around QTP show positive anomalies (Figure 3), which are in agreement with the statistics between the Bouguer gravity anomaly and the topography (Woollard, 1959). According to the Airy isostatic hypothesis, the crustal root is necessary under the huge mountains, and the relatively low density of the crustal material instead of the high-density mantle material under the mountain may lead to mass deficiency, which shows obvious low anomalies on the Bouguer gravity map. Therefore, under normal circumstances, the high elevation shows low anomalies, and the low elevation corresponds to high anomalies for the Bouguer gravity. For this reason, the SB as the biggest basin in this area shows significant Bouguer positive gravity anomalies because of the relatively upward bulge of the Moho surface to the QTP .

The Top of the Crust
The average depth of the field source of 1-D is approximately 5 km (Figure 4), which is nearly the top of the upper crust. On the 1-D map, there are several areas with intense gravity changes on the margin of QTP, including QLB, LSB, NSSB, part of QTB, LB, and HB. In these areas, the high and low gravity anomaly belts are arranged alternately, and the directions of these anomaly belts are highly similar to the local topographic relief. To intuitively compare the consistency between the terrain and the anomaly belts, we used the hydrology tools in ArcGIS (Jenson & Domingue, 1988) to extract ravine lines based on DEM data. The similarity between the directions of the valleys and anomaly belts (yellow boxes in Figure 4) indicates that the development of ravines in this area is mainly controlled by the geological structure at the top of the upper crust.
At this depth, low-gravity anomalies mean the local low-density anomalies on the top of crustal, which might be caused by low-density material such as the basin sediment or the fault fracture zone. However, the gravity anomaly belt at this depth is all narrow and long shaped, which is more likely the reason for the fault fracture.
In the LSB region, the major axis directions of positive and negative gravity anomaly belts are NW and NE. The NE anomaly belts are parallel to the Longmenshan fault (F6), and the NW anomaly belts are the structures developed during the Indo-Chinese orogenic phase and the Yanshanian orogenic phase (Tian et al., 2017).
The north subblock of the SYB, the NSSB, also shows obvious gravity anomaly belts on the 1-D map, which are mainly distributed in the SN direction and are bounded by the Lijiang Xiaojinhe fault (F11). On the south

10.1029/2019EA000943
Earth and Space Science of F11, the gravity anomalies change gently without obvious long-axis direction in the MYSB area. The huge difference between the two subblocks of the SYB may imply that it is not reasonable to regard the SYB as a whole. The gravity anomaly belts in QLB are several NW strikes parallel to the North Qilian Fault (F3) and are very straight. To the west of the Jinshajiang fault (19), which is distributed along the valley of the Jinsha River, there are two remarkable low-gravity anomaly belts that extend along two other large rivers of the Sanjiang area, the Lancang River and the Nujiang River. The strikes of these belts gradually change from the NW direction at the northwest end to the nearly the SN direction at the south end.

The Upper Crust
The average field source depth of 2-D is approximately 9 km ( Figure 5), which is about the upper crust. On the 2-D map, the frequency of abnormal signals decreases, and the alternating positive and negative gravity anomaly belts disappear. At the edge of the QTP, the intensity of the gravity anomalies is still greater than in other regions. In the OB, SB, and AB blocks, the gravity changes are much more gentle than in other tectonic regions. On the northwest side of OB, there is a significant low-gravity anomaly belt between two of the Peripheral Ordos faults (F20). This anomaly consists of a NE trending part on the west and a near-EW trending part on the east, which corresponds to the Yinchuan-Hetao basin by topography.

The Middle Crust
The average field source depth of 3-D is approximately 25 km ( Figure 6), which is about the middle crust. On the 3-D map, the scale of the gravity anomaly increases, and the frequency decreases. The mountain areas, such as QLB and LSB, mainly show obvious low-gravity anomalies, which are opposite to the topography. However, the topography and the gravity anomalies show synchronous changes along the Yinchuan-Hetao basin in the NW corner of the OB. The Yinchuan-Hetao basin is low in both topography and gravity anomalies.
At this depth, the low-gravity anomalies on the detail map indicate a material density deficit, which may be due to crustal flow in this area . The flow capability of the crustal material may be related to the increase in fluid, and the increase in fluid material corresponds to the anomalously low density of the material. Yang et al. (2015) compared the magnetotelluric resistivity profile (Jin et al., 2010) with density anomalies at corresponding depths in the middle of the QTP and confirmed that the low-density anomalies were coincident with the low-resistivity anomalies. Therefore, based on density anomalies derived from gravity data, Yang et al. (2015) discussed crustal flow locations, including bottom splitting upwelling locations and lateral flow locations. The low-density zone on the southern margin of the QTP is very wide and may be the source of the lateral extrusion of matter (red triangles in Figure 6) , while the low-density anomalies (low-gravity anomaly area) located in west Qinling and the SYB around the SB represent the channel locations of crustal flow. Although it is impossible to discuss the time when the flow itself started, gravity anomalies can reflect the locations and scales of these channel flows.
The lower crust flow explains the difference between gentle and sharp changes from the Qinghai-Tibet Plateau to the surrounding topography (Clark & Royden, 2000;Royden et al., 1997), but the topography still rises within the range affected by the crust flow. This theory cannot explain the corresponding relationship between low topography and low density in the Yinchuan-Hetao Basin on the northeastern margin of Ordos. Wang et al. (2017) believes that the preexisting high-density crustal root in this area is the reason for the obvious large scale of the low-gravity anomaly in the basin without topographic uplift.
In SYB, the south and north subblocks show significant differences on the 3-D map. NSSB mainly consists of low-gravity anomalies, while MYSB mainly consists of high-gravity anomalies, and the boundary of the anomaly division coincides with the Xiaojinhe-Lijiang fault (F11). Around the MYSB, the low-gravity anomalies are coincident with the two low crustal channels in this area (Bai et al., 2010;Yang et al., 2015;Li, Liu, et al., 2017;Fang & Hou, 2019), as indicated by the two arrows in the south on Figure 6. The high-gravity anomalies in the MSYB and the channel flow around it, indicating that this subblock blocks the flow of the lower crust, just like the SB does but with less solidity. The Qinling Mountain (along the  Figure 6. According to Yang et al. (2015), the crustal flow channel in the north distributes along the west Qingling and the Yinchuan basin, as indicated by the northeast pointing arrow in Figure 6. However, in the result of the 3-D map (Figure 6), the crustal flow toward the north is blocked by some high-gravity anomalies in this channel.

The Lower Crust
The average field source depth of 4-D is approximately 49 km (Figure 7). On the 4-D map, the contours of different tectonic units are clearly reflected and are limited by major faults. The earthquakes in the area are also concentrated in the transition zones of different tectonic blocks. SB is characterized by obvious high-gravity anomalies, while the Longmenshan secondary block on the west side is a low-gravity anomaly area. NSSB in the northern part of SYB is a low-gravity anomaly area, while MYSB in the southern part is a high-gravity anomaly area.
Along the Xianshuihe fault (F7), the isolines of the low-gravity anomalies of LSB and NSSB show homomorphic distortion, meaning that this fault is large and deep enough to affect the crust density structure in this area, but the structures on both sides are similar in density. Furthermore, the crustal structures reflected by the 4-D map of gravity anomalies are very similar along the Longmenshan fault (F6) and the Xiaojinhe-Lijiang fault (F11), with a low-gravity anomaly in the Northwest and a high-gravity anomaly in the Southeast. The gravity anomaly contrasts on both sides of these faults may be due to the differences in crust thickness, which are derived from the blocking of the eastward flow of QTP material by hard blocks such as SB and MYSB. The differences in crust structure on the two sides of this boundary can also be observed in the crustal velocity structure model , meaning that the Xiaojinhe-Lijiang fault is very similar to the Longmenshan fault in the geology settings and should be noticed in the seismic hazard research in this region.
The Moho depth in the QTP region is approximately 60 km (Bao et al., 2015;Shin et al., 2015;Lu et al., 2019), which is close to the average source depth on the 4-D map (approximately 49 km). Therefore, the anomalies on the 4-D map reflect large-scale changes in the crustal structure, both in the crustal thickness and in the crustal material density. The gradient zones of gravity anomalies are often considered the boundaries between tectonic units and provide important tectonic background for strong earthquakes. Although most earthquakes in this area are distributed in the upper crust, the occurrence of earthquakes is still affected by the deep tectonic background, and the clustering phenomenon of earthquakes in the gravity anomaly conversion zone also strongly supports this point. For these reasons, the 4-D map is chosen as the basis for the study of the numerical relationship between gravity changes and earthquake frequencies in this research.

Discussion
The fourth-order wavelet detail map (4-D) can clearly indicate the outlines of the basic tectonic units in this area. The tectonic units located in the margin of the QTP, such as the LB and HB in the south, NSSB in the southeast, LSB in the east, and QLB in the north, show obvious low gravity anomalies, while the tectonic units located in the periphery of the QTP mainly show high-gravity anomalies, such as the SB in the east and TRB in the north. The gravity anomalies corresponding to these tectonic units can roughly reflect their positions and ranges, indicating the differences between tectonic units in terms of the crustal structure ( Figure 7). Fault belt tectonic units such as the Altun fault belt (F1), the north Qilian fault belt (F2), and the Longmenshan fault belt (F6), which belong to the northern and eastern boundary tectonic zones of the QTP (Deng et al., 2002;Deng et al., 2003), show the transition zones between high-and low-gravity anomalies, which are also the locations of earthquake clusters (Figure 7).
To accurately discuss the degree of this spatial clustering, we need to quantitatively measure the spatial variation of gravity anomalies on the 4-D map and the statistical relationships between earthquake distribution and gravity anomaly variation. Finally, in this section, a simulation of seismic risk assessment is preceded by dividing the earthquake catalog into two parts at 2008 to test the effectiveness of the assessment using this distribution relationship.

Quantitative Description of the Spatial Variation of Gravity Anomalies
The spatial variation characteristics of the gravity field are usually described by gradients in three directions, that is, x, y, and diagonal directions, and this field is only sensitive to the tectonic structure of a particular direction. In this paper, we need to construct indicators that can describe the spatial variety of gravity in any direction as a proper statistic for further research. Based on the 4-D map, we calculated four different gravity parameters as statistically independent variables based on the slope and the curvature of the gravity anomaly values.
Slope is the most commonly used topography index to characterize changes in the spatial terrain. Regardless of the unit and physical meaning, by substituting gravity for elevation, we obtain the gravity slope factor, which is the maximum rate of change in the gravity value from one cell to its neighbors. For each grid cell on the 4-D map, the slope is (Burrough & Mcdonell, 1998): For the center cell e in Figure 8, the 3 × 3 cell neighborhood is used to calculate the change rates in the x and y directions: where x_cellsize is the cell size in the x direction and y_cellsize is the cell size in the y direction.
On the SLP map, some earthquakes are clustered on the edges of large-slope zones, especially the 2001 Ms8.1 Kunlun Earthquake (Figure 9). To describe these positions, three more factors related to curvature are also

10.1029/2019EA000943
Earth and Space Science constructed. For each grid point, the curvature K p of profile p in an arbitrary direction can be calculated as follows: The gravity curvatures in 12 directions with a 15°angle between two adjacent directions are written as K p1 , K p2 ,⋯,K p12 , and the factors of the gravity curvatures are written as follows: These gravity spatial change factors represent the tectonic backgrounds of earthquakes, and the Bouguer gravity anomaly data used to calculate these factors come from a static model, the GECO model, which illustrates the static part of the gravity field related to the existing crustal structure. The structure of the crust gradually changes with earthquake activity. Tanaka et al. (2001) noted that the coseismic gravity change for a M6.1 earthquake was 6 μGal. shown in the maps were between −100 and 60 and between −140 and 140 μGal, respectively. However, the Bouguer gravity anomaly reflecting the crustal structure changes in the study area was between −700 and 30 mGal. According to this estimation, the Bouguer gravity anomaly change caused by earthquakes during one century may not be large enough to cause significant variation in the gravity spatial change factors, including the SLP value. Therefore, these numerical factors obtained based on Bouguer gravity anomalies also have sufficient stability, and their changes over decades can be ignored when using them to study the seismic tectonic background.

Quantitative Statistics of Gravity Anomaly Changes and Earthquake Distributions
The earthquake frequency distributions of different levels of the four factors were used to study the relationship between the earthquake distribution and the gravity field characteristics. For each factor, the value was divided into 100 or 20 classes. To remove the impact of the area of different value ranges, the quantile method was used during the reclassification of attribute maps. The quantile method made each class contain a similar number of cells. Then, the cell numbers of each value range N i c were used as the denominators of the frequency of earthquakes to completely eliminate the impact area. Finally, the standardization of the earthquake frequency was used to compare the different factors. The parameter f n (i) is the earthquake frequency of the i class of any factor, F n (i) is the normalized f n (i), and N i e is the earthquake number within the i class area on this factor map. The formula can be written as follows: The bivariate correlations of the factor classes with the earthquake frequency F n are calculated using the Spearman correlation coefficient (Tables 1-3). The correlation coefficients of different earthquake subsets and different factors are shown in Figure 13.
To answer the questions of whether the statistical characteristics can be affected by the reclassified number of factors, the aftershock removal, and the appearance of large earthquakes for which the recurrence period is much longer than the earthquake catalog, we divide the earthquake catalog into 12 subsets to analyze the influences of these conditions. Subsets 1 to 6 form the group with 100 classes after reclassification, while subsets 7-12 form the group with 20 classes. In each group, the first three subsets, that is, 1-3 and 7-9, are the teams that are all recoded before aftershock removal, while the last three subsets, that is, 4-6 and 10-12, are the teams with the main shock determined after aftershock removal. In each team, the first subset, that is, 1, 4, 7, and 10, is the earthquake catalog during the whole time period from January 1970 to June 2019; the second subset, that is, 2, 5, 8, and 11, is the earthquake catalog before 2008 ( 1970.01-2007.12); and the third subset, that is, 3, 6, 9, and 12, is the earthquake catalog after 2008 (2008.01-2019.06).
Based on the above grouping, the earthquake catalogs of magnitude above 3.0, 4.0, and 5.0 are studied to analyze the influence of the magnitude on the correlation coefficients between the earthquake frequency and gravity factors (corresponding to the solid line, dashed line, and dashed dotted line in Figure 13, respectively). As the sample size of the catalog of earthquakes with magnitudes of 5.0 or above is too small, only the 20-classification group is counted (the dashed dotted lines in Figure 13 cover numbers only 7-12).
The correlation coefficients of the 20-classification group are significantly higher than those of the 100-classification group. The reason behind this difference may be the sample size, which reflects the average number of earthquakes in each grade (as shown by the gray lines in Figure 13). When the average sample size for each grade is greater than 200, such as the team with all records in the 20-classification group (red solid line between 7 and 8 in Figure 13), changing the sample size no longer causes a significant change in the correlation coefficient. Approximately half of the correlations in the subset of earthquakes above M5.0 are not significant. This feature may be due to the rapid decrease in the number of earthquakes above M5.0, for which the average sample size of each level is also less than 50. The small sample size makes it difficult to reflect the real statistical characteristics. Therefore, a sufficient average sample number of each grade may be necessary to obtain a more reliable correlation analysis.
The period of the earthquake catalog also has a significant impact on the correlation coefficient, especially in the teams assessed before aftershock removal (i.e., 1-3 and 7-8). The lines of the three curvature factors (AVE, MAX, and MID) show significant sharp downward corners for these teams (blue, purple, and green lines, respectively, between 1-3 and 7-8 in Figure 13). Although there are more earthquake catalogs before 2008, the subsets after 2008 have higher correlations. This phenomenon is not obvious for the SLP factor, except for the subset of teams without aftershock removal and with insufficient sample sizes (the red dashed line between 1 and 3 and the red dashed dotted line between groups 7 and 8 in Figure 13). This change may be related to the Wenchuan earthquake in 2008, and after removing aftershocks, this phenomenon can be well suppressed.
As the minimum magnitude of the earthquakes in the catalog increases, the correlation coefficient seems to decrease. However, the sample size decreases rapidly as the minimum magnitude increases, which also leads to a decrease in the correlation coefficient. Upon discarding the statistical results for M5.0, in which the sample size is too small to yield a reliable correlation coefficient, by comparing the results of SLP in the subsets of M3.0 and M4.0 (solid and dashed red lines in Figure 13, respectively), the teams of main shocks show that an increase in magnitude leads to an increase in the correlation (between 10 and 12 in Figure 13), while the teams without aftershock removal show the opposite trend (between 7 and 9 in Figure 13). In connection with the suppression effect of aftershock removal on the time impact, aftershock removal is clearly necessary in the statistical analysis of correlations between gravity change factors and earthquake frequencies.
Based on the above description, a step plot was drawn for the subset of the main shock team in the 20-classification group (subset number 10 in Table 1) to observe the relationships between different attribute factors and the distribution of the earthquake frequency. Only the SLP factor was used in the initial magnitude comparison (the SLP of subset number 10 in Tables 2 and 3, which has the highest correlation coefficient). According to the reclassification value, the earthquake frequency was used to analyze the earthquake distributions at different attribute levels. The step plots of the earthquake frequency relative to different factors are shown in Figure 14. In the 4-D step plot, two value ranges for high earthquake frequency are obvious. The gravity values of these two ranges are −17-−1 and 13-15 mGal, while the gravity values of the whole 4-D map vary from −42 to 55 mGal. The step plot of MAX shows that the earthquake frequency increased with the maximum of the directional curvature, but the Spearman correlation coefficient was only 0.563, which was significant at the 0.01 level. This value is not as good as those of the other two curvature factors, AVE and MID. The best directional curvature factor is for the AVE, with the highest correlation coefficient of 0.716, which is significant at the 0.01 level.
For the SLP, the earthquake frequency shows an obvious increase with SLP from 11°to 18°and otherwise decreases through the whole range of 0°-30°. The Spearman correlation shows that the earthquake frequency increases with the SLP, which has coefficients of 0.878, 0.933, and 0.764 for magnitudes above M3.0, M4.0, and M5.0, respectively, with significance at the 0.01 level based on a two-tailed test. The correlation coefficients increase from 0.878 to 0.933 when the initial earthquake magnitudes increase from M3.0 to M4.0, meaning that stronger earthquakes are prone to clustering in steeper-slope areas, where the gravity changes sharply. The decrease in the coefficient to 0.764 for the subset with initial earthquake magnitudes of M5.0 may be due to the decrease in the sample size.

The Efficiency of the SLP Factor in Earthquake Risk Assessments
The Wenchuan Ms8.0 earthquake was the largest earthquake in this region in decades and has had great impacts on society and the economy. At the same time, the location of the earthquake is also a place where the SLP values are relatively large (17.8°), which is consistent with our hypothesis that the gravity gradient area should be a region with a high earthquake risk because of the high SLP values, which provide a digital expression of the spatial changes in crustal structure. However, the Longmenshan Fault was seismically quiet before the Wenchuan earthquake, which resulted in an underestimation of the earthquake risk (Xu et al., 2009). In this case, can we determine the correlations between the SLP parameters and the earthquake frequency? The answer is yes. The earthquake subset before 2008 showed a significant positive correlation between the SLP class and earthquake frequency (the correlation coefficient of the subset after the aftershock removal is 0.848; see the SLP column in row 11 of Table 1). However, to indicate the effectiveness of SLP parameters in earthquake risk assessments, further discussion requires statistical methods. Therefore, for the SLP with a high correlation coefficient with the earthquake frequency before 2008, we carry out a seismic risk assessment based on SLP values in the eastern QTP. Taking 2008.01 as the starting time of the declustered earthquake catalog, the effectiveness of the SLP in earthquake risk assessments is tested in this section.
To quantitatively evaluate the seismic spatial assessment capacity of the gravity factor SHP, receiver operating characteristic (ROC) curves are drawn for the earthquake catalog above M3.0, M4.0, and M5.0. Taking each cell in the raster data set as a prediction sample, if any earthquakes in the earthquake catalogue are located in the cell, then this cell is considered a positive event; otherwise, it is recorded as a negative event.
For each class i of SLP, the cell with the SLP class larger than i is considered as a positive prediction event; otherwise, it is predicted as a negative event.
For each SLP class i, the cells with earthquakes inside and SLP classes larger than i are true positive events (TP), those with earthquakes inside and SLP classes no larger than i are false negative events (FN), those without earthquakes inside and SLP classes larger than i are false positive events (FP), and those without earthquakes inside and SLP classes no larger than i are true negative events (TN). The false positive rate (FPR) and the true positive rate (TPR) can be written as shown below: Then, the TPR and the FPR are calculated for each class of SLP. Taking the FPR as the x axis and the TPR as the y axis, the ROC curve can indicate the assessment efficiency of the SLP for earthquake risk spatial assessments ( Figure 15).
The more the ROC curve deviates from a diagonal line, the better the assessment efficiency is. According to the ROC curve, the assessment efficiency increases with the initial magnitude of the earthquake catalog. To

10.1029/2019EA000943
Earth and Space Science numerically evaluate the assessment efficiency, the area under the curve (AUC) of each ROC curve is calculated ( Figure 15). The AUCs increase significantly from 0.580 for the earthquake catalogue above M3.0 to 0.614 for that above M5.0. However, earthquakes decrease exponentially as the magnitude increases, and the statistical stability of the ROC shows an obvious deterioration for the earthquake catalogue above M5.0. A comparison of the three ROC curves and their AUCs shows that stronger earthquakes correspond to higher assessment efficiencies of SLP, which means that strong earthquakes tend to occur in regions with sharp spatial changes in gravity anomalies. The results in this study support the view that strong earthquakes tend to cluster in regions with a spatial change in the crustal material density, which can be derived from the low-frequency component of the gravity field.
Although the ROC curve shows that the phenomenon of earthquake clustering along the gravity gradient zone is real, some problems remain in SLP-based assessment according to the SLP map ( Figure 9). On the one hand, some earthquake clusters are underestimated. A number of earthquakes smaller than M5.0 are located in part of WYB and QTB, where the SLP class is less than 10. This type of failure may be because the gravity anomalies in 4-D only occur around active faults with enormous changes in the crustal mass on both sides in large scales. For this reason, the SLP map derived from 4-D may miss areas with young and small active faults without enough crust mass differentiation.  On the other hand, there are areas with high SLP values but without many earthquakes, such as the high-class belt along the Eastern Kunlun Fault (EKF, F14 in Figure 9) and the Qinling orogenic belt (south of F21 in Figure 9). For the EKF belt, this phenomenon may be due to the stress release during the 2001 Ms8.1 Kunlun Earthquake. For the Qinling orogenic belt, it may be related to the similar directions of block movement on different sides of the Qingling, which cause only slight relative movements along this belt (Zhao et al., 2012).

SLP Distribution and the Spatial Variation in Crustal Structure Across the Eastern QTP
On the eastern boundary of the QTP, two high SLP value belts are very striking along the Longmenshan fault belt (F6) and the Xiaojinhe-Lijiang fault belt (F11). In these zones, the gravity changes rapidly across the fault belts, which may indicate rapid changes in the crustal structure, thickness, or density. On the south side of the Xiaojinhe-Lijiang fault belt (F11), there is an obvious high-gravity anomaly center originating from the high-density anomaly of the MYSB block, which has some similarity with the SB. Correspondingly, the region is also characterized by a high Love wave group velocity (Zheng & Wang, 2017), high Rayleigh wave velocity (Bao et al., 2015), high P-wave velocity , and high resistivity (Cheng et al., 2015). These high-density and high-velocity anomalies indicate a stronger crust under the MYSB block, which blocks the low-density and low-velocity weak material of the crustal flow and separates the flow into two branches (Figures 6 and  7). The situation around the MYSB block is very similar to that around the SB block, which also contains significant high-velocity anomalies and separates the low-speed crustal flow material into two branches (Bao et al., 2015). The similarity of the high-gravity anomalies of these two blocks and the low-gravity anomalies of the crustal flows around them are also very obvious on the 4-D map (Figure 7). Zheng and Wang (2017) speculated that the intrusion of mantle-derived materials caused the abnormal density and mechanical strength in this area, which were greater than those of the surrounding area and caused the blocking. The structural characteristics of this area can be compared with those of the SB block and the Longmenshan fault belt (F6). The intersection between the high-strength blocks on the east side and the crustal flow on the west side resulted in rapid uplift and steep relief, which also mean this area has a high risk of strong earthquakes along the Xiaojinhe-Lijiang fault belt (F11).

Conclusions
To study the changes in crustal structure based on gravity anomalies and their relationships with the spatial distribution of earthquake frequency, we analyzed the horizontal changes of crustal structure based on gravity data and derived four gravity factors, SLP, MAX, AVE, and MID, from the 4-D map to indicate the spatial gravity field changes.
The results of the statistical analysis show that the SLP, AVE, MID, and MAX show significant positive correlations with the earthquake frequency at the 0.01 level and have relatively high correlation coefficients of 0.878 (SLP), 0.716 (AVE), 0.710 (MID), and 0.563 (MAX) (for the declustered subset with 20 classes from 1970-2019 above M3.0, i.e., number 10 in Table 1), meaning that greater variations of the gravity field correspond to Figure 14.
Step plots of earthquake frequency in different factor classes. The horizontal axis is the classification number of factors, and the vertical axis is the normalized value of the earthquake frequency in each factor class.

10.1029/2019EA000943
Earth and Space Science a higher frequency of earthquakes, and earthquakes tend to occur in areas with spatial changes in the crustal gravity anomaly gradient.
Based on SLP, the index with the highest correlation coefficient, a test of the spatial assessment of earthquake risk is performed. The ROC curve shows that the SLP is useful for assessing the spatial earthquake risk. Contrasting the AUCs of the ROC curves of the earthquake catalogue above M3.0 (0.571), M4.0 (0.598), and M5.0 (0.627) clearly shows that the assessment efficiency increases with the initial magnitude of the earthquake catalog. The results indicate that the stronger the earthquakes are, the more likely the earthquakes are to cluster in regions with spatial changes in gravity anomalies on a 4-D map.