Surface Slip Distribution and Earthquake Rupture Model of the Fuyun Fault, China, Based on High-Resolution Topographic Data

The characteristics of earthquake surface ruptures, such as geometry, slip distribution, and coseismic deformation, contain important information about the earthquake rupture process, and so investigations and analyses of earthquake surface rupture have played a crucial role in modern earthquake hazard studies, especially with the increasing availability of high-resolution topographic and imagery data for tectono-geomorphic interpretation. In this study, we use Structure from Motion (SfM) photogrammetry to build a 1m resolution digital elevation model (DEM) of the fault and combine this with filed observations to map the surface ruptures of the 1931 M8.0 Fuyun earthquake, China. These high-resolution topographic data enable to identify and measure the displaced gullies, and so the rupture locations and along-strike slip distribution are obtained in detail. Four paleoearthquake events are identified through the offset cluster characteristics. The coseismic offset of the 1931 Fuyun earthquake is extracted from the offset distribution, which shows four continuous undulations along the fault strike, corresponding to the four segments of surface rupture. Moreover, a high offset gradient is observed in the step area connected by the rupture segment. These findings, combined with the width and bending angle of the step area at the joint of the rupture segment, indicate that the 1931 Fuyun earthquake was a cascade rupture formed by four rupture segments. This study expands the available offset measurement data of Fuyun fault and confirms the applicability of high-resolution topographic data to active tectonic research.


Introduction
Moderate-to-large earthquake (Mw 6-7) could rupture the ground surface and produce various displacement landforms which can be used to constrain along-fault slip distribution [1]. The characteristics of earthquake surface ruptures, such as geometry, slip distribution, and coseismic deformation, contain important information about the earthquake rupture process, and so investigations and analyses of earthquake surface rupture have played a crucial role in modern earthquake hazard studies [2]. By investigating the displacements along the fault zone, we can better understand critical aspects of the earthquake process, including source complexity [3][4][5], rupture propagation dynamics, and slip histories over multiple earthquake cycles [6][7][8]. Moreover, such information can also be used to predict the future behavior of the fault, which is beneficial for seismic risk assessment and loss mitigation [9]. Many early studies were based on traditional field mapping as well as detailed local topographic surveys and air photo interpretation [10]. However, this method was time-consuming, and the subsequent analysis of earthquake faulting is also restricted to scarce data [11]. High-resolution topographic data have now become increasingly available, and their applications have significantly improved geologist's abilities to map and measure surface rupture and deformation. For example, airborne LiDAR has been successfully used to map the rupture and vertical slip and provides a new for topography mapping with high resolution [12,13]. Compared with traditional field survey methods, LiDAR enables to take surface faulting sampling at much higher spatial densities than field surveying, showing the advantage of LiDAR to survey the extensive and complex rupture associated with the earthquakes [14]. However, the high cost in equipment and data acquisition has greatly restricted LiDAR applying in geoscience programs [15].
In contrast, the low-cost and flexible aerial photogrammetry known as the Structure from Motion (SfM) technique is suitable for acquisition of topography rapidly in large areas [8,16]. Unmanned aerial vehicles (UAVs) are recently gaining increasing popularity in aerial photogrammetry methods, since the autopilot technology to flight control and image acquisition. Furthermore, UAVs are able to fly at a much lower altitude and thus can collect images with an unprecedented level of detail and provide a great opportunity for high-resolution mapping, which not only improve the identification and measurement of geomorphic offset, but also contribute to understand critical aspects of the earthquake process. It has now been employed in a number of structural geomorphology applications, which have obtained new insights in fault activity research and earthquake process analysis [5,11,17,18].
We explore the application of SfM photogrammetry for identifying surface ruptures and establishing fault offset distributions along the Fuyun fault on which the 1931 M8.0 Fuyun earthquake occurred. The Fuyun fault is a large right-lateral strike-slip fault zone in the Altay Mountains in Northwest China. The 1931 Fuyun earthquake led to large-scale surface rupture of 159 km and affected a wide area [19]; as one of only a few very large earthquakes in the world at that time, it attracted substantial attention from the seismological community. Due to the sparsely populated area, minimal human activities, and favorable natural conditions, various geomorphic phenomena are still well preserved almost a century later, making this region a natural laboratory for investigating surface rupture using SfM photogrammetry. Klinger et al. used Quickbird satellite to obtain a panchromatic image of the Fuyun fault with a resolution of 0.6 m [20]. In this study, we provide new high-resolution topographic data (1 m) along the Fuyun fault obtained by UAV. The geomorphic indicators of offset gullies along different fault segments are identified and carefully measured. Based on statistical analysis of displacement measurements, we cautiously restore the history of strong earthquake rupture and propagation characteristics along the fault. Our data are critical for informing paleoseismic, tectonic geomorphology, and structural geologic mapping of the Fuyun fault, as well as for investigating the behavior of the Fuyun fault by limiting the slip distribution along the fault trace.

Geologic Setting
2.1. Regional Tectonic Environment. The Altay orogenic belt is an active seismic belt located in the interior of Central Asia (Figure 1), at the geographic confluence of China, Kazakhstan, Russia, and Mongolia and the tectonic location between the north-south compression and convergence area of Tien Shan mountains and the Baikal NW-SE extensional rift. Although this region is far from the collision zone between the Indian plate and Eurasian plate, it remains an inland active zone within the Eurasian plate.
In the late Paleozoic, due to amalgamation and compression of the Mongolia-Siberia plate and the Junggar-Kazakh plate, a strong orogeny occurred in the splicing zone, which formed the Altay orogenic belt. This mountainous belt is composed of a series of arc-shear-nappe tectonic systems [22]. During the Cenozoic, the tectonic activity of the Altay orogenic belt began to strengthen and was mainly characterized by fault-block uplift and differential movement, which activated early ductile faults and led to brittle deformation, as evidenced by the Irtysh fault and Mayin-Erbuo fault, and reformed a series of NNW trending faults, including the Fuyun fault. During the Quaternary period, tectonic movement was further strengthened, increasing the scale of the Altay orogenic belt [23].
The structure of the Altay Mountains is affected by both the remote effect of collisional deformation of the Indo-Eurasian plate [24] and the dynamics of the Mongolia-Siberian plate. This enhances regional activity, which is characterized by a large right-lateral shear fault with a NNW strike and a right-lateral thrust fault with a NWW strike [25]. Many earthquakes have occurred in the Altay region over the past hundred years: a M8.0 earthquake occurred in Fuyun, Xinjiang, China, in 1931 [26], a Mw8.1 earthquake occurred in Altay, Gobi, Mongolia, in 1957 [27], and a Ms7.3 earthquake occurred in Zaysan, Kazakhstan, in 1990 [28]. In addition, several small earthquakes have occurred in this region. However, in front of the Altay Mountains, only intermountain basins controlled by strikeslip faults have developed; these basins are small in scale, unlike the large foreland basins in front of the Tien Shan mountains or on the periphery of the Qinghai-Tibet Plateau.

Fuyun
Fault. The Fuyun fault, also known as the Koktokay-Ertai fault, is located at the southwest foot of Altay Mountains and the northeast margin of Junggar Basin. The strike is generally NNW, starting from the Maizengsayi Mountains in the north and extending southward along the southwest edge of the Altay Mountains to the banks of the Ulungur River ( Figure 2). The administrative divisions including this fault are Fuyun County and Qinghe County, Altay City, and the Xinjiang Uygur Autonomous region. The Fuyun fault zone trace includes several faults with NNW and almost NS strikes. From north to south, the fault strikes almost north-south and then bends slightly to the southwest, before extruding slightly eastward, exhibiting a gentle wave in the planar dimension, with an overall strike of 342°, a dip angle of 60-85°to the NE, a total length of almost 200 km, and a width of tens of meters to more than 100 m [29].
Since the neotectonic period, the southern foothills of the Altay Mountains have predominantly exhibited horizontal movement and developed asymmetric elbow-shaped river systems, which are also clearly reflected in the Fuyun fault. Fault landforms such as offsetted river systems, mountain ridges, inverted rock piles, and fault scarps are particularly well developed along the fault. Obvious dextral    The Fuyun fault area has a typical temperate continental climate, with long and cold winters and short and dry summers. The average temperature from 1981 to 2010 was 3.8°C, and the average annual precipitation was 209.7 mm (China Meteorological Data Network, https://data.cma.cn/). The dry climate effectively preserves geomorphic features such as surface ruptures and gullies. In addition, the study area is sparsely populated and minimally affected by human activities, which is conducive to studying the geomorphology of this region.
2.3. Background to the 1931 Fuyun Earthquake. The Fuyun fault zone is located in the interior of the Asian continent, in the Pamir-Baikal seismicity area of western Mongolia. Since the last century, this area has become one of the most seismically severe intracontinental regions in the world [30]. According to incomplete historical statistics of the past 400 years, there are more than 420 earthquakes with M ≥ 5:5. Especially since 1900, earthquakes with M ≥ 8 have occurred seven times [29]. At approximately 05 : 00 on August 11, 1931, an earthquake of M8.0 occurred in Fuyun. The macroscopic epicenter was located in Karaxingar (46°44.5′N, 89°54.0 ′ E) [29]. The seismic waves recorded by the San Juan Seismological Station in Puerto Rico lasted for two and a half hours. The earthquake affected a wide area covering thousands of kilometers from Semipalatinsk in Russia in the west to Lanzhou in the east. The earthquake was so strong that vast areas of China and Mongolia, including Urumqi, Haba River, Balikun, and Gilgrantu of Mongolia, suffered varying degrees of damage.
Surface damage from this earthquake was predominantly limited to the narrow zone between Koktokay in Fuyun County and Qinghe County and characterized by a series of fault scarps, grooves, and ridges, accompanied by collapses and landslides [26]. A surface rupture zone formed 159 km along Koktokay-Ertai fault with NNW strike [19]. According to the intensity distribution of the Fuyun earthquake in China, combined with overseas intensity data, Shi et al. [31] stitched together a complete image of the earthquake crack distribution (  4 Lithosphere reached XI on the China seismic intensity scale. The east side of Karaxingar was seriously damaged, forming a spectacular subsidence area with vertical and horizontal internal cracks, generating a seismic fault zone passing through the epicenter area, splitting mountain ridges, cutting off river systems, forming dammed lakes, and causing the collapse of rock piles and vegetation. Before the Fuyun earthquake, there were obvious foreshocks in the earthquake area, and the seismicity increased significantly between 1917 and 1927, with eight successive moderately strong earthquakes. After the earthquake, aftershocks were frequent and lasted for a long time. The largest aftershock occurred on the seventh day after the main shock with a magnitude of 7 1/4. The epicenter was located near Tosbastao, approximately 40 km south of the epicenter of the main shock. The other major aftershocks predominantly occurred in the southern part of the fault zone [29].
Regarding the seismogenic structure of the 1931 Fuyun earthquake, the Seismological Bureau of the Xinjiang Uygur Autonomous region suggested that it occurred in the Fuyun fault zone [29]. On the one hand, the seismic fault zones produced by this earthquake typically exhibit rightstepping en-echelon formations, feathered structures, and cross-cut ridges and gullies, forming a horizontal fault distance of several meters, reflecting the characteristics of right-lateral horizontal deformation. On the other hand, the region corresponding to a value of X on the China seismic intensity scale is distributed along the Fuyun fault in a narrow band, and the long axis of the isoseismal line is consistent with the strike of the fault. The seismic fault zone is also limited to the Fuyun fault and its southward extension area, with no seismic fault traces found in the nearby Ertix fault. Both new and old faults exhibit the same type of movement, indicating that the Fuyun fault was the seismogenic structure of the 1931 earthquake, and that this large earthquake was a manifestation of new activity along the fault, as well as there is a new development to the south [26].
In the 1980s, field investigations measured strike-slip offsets of 6-10 m [29]. In recent years, Klinger et al. [20] obtained the offset distribution of the Fuyun fault by interpreting remote sensing images, which revealed a coseismic offset in 1931 of 6:3 ± 1:2 m, as well as multiple other offsets related to the 1931 event, indicating at least four distinct previous earthquakes. These slip distributions were similar for all five consecutive earthquakes, indicating that the Fuyun fault conforms to the characteristic seismic model. However, remote sensing images are two-dimensional; therefore, fluctuations of geomorphology cannot be observed. Moreover, due to the limited data quality, the recognition of broken landforms is greatly limited. Finally, human subjectivity brings substantial uncertainty to the displacement measurement; thus, further analysis of the rupture process and seismic risk is hindered.

Methodology and Data
3.1. Principle of SfM Photogrammetry. The SfM method was originally used in the field of computer vision [32]. By tracking and photographing moving objects, the three-dimensional structure and motion of moving objects can be restored from multiview two-dimensional images. With the development of computer vision theory and efficient automatic feature matching algorithms, SfM technology has become widely used in digital photogrammetry for the three-dimensional entity reconstruction of highly overlapping images [33]. Digital photogrammetry produces stereoimagery by taking two pictures of the same object from different angles. According to the geometric inversion theory of the photography process, two-dimensional points are projected and transformed to reconstruct the three-dimensional model [34,35]. This measurement method requires the known camera position and motion to determine the three-dimensional coordinates of the image. Moreover, the SfM method is based on the principle of Scale-Invariant Feature Transform (SIFT), which obtains images with sufficient overlap from multiple angles, and performs automatic matching of image features, which greatly improves the degree of photogrammetry automation [36]. By correcting the bundle adjustment of matching features, the XYZ coordinates of each feature in the image can be identified, and the three-dimensional position and direction of the camera are automatically solved to form a sparse point cloud frame [37]. In order to further improve the point cloud density, searching and matching between images is conducted based on the principle of multiview stereopsis (MVS) to obtain dense point cloud data [38].
Because the 3D terrain model generated by the SfM method does not contain spatial location or scale information, it is necessary to correct its geographic location and scale. Traditional georeferencing based on ground control points (GCPs) provides reliable positioning; however, the geometric accuracy depends on the number and spatial layout of GCPs, which are limited by time and cost constraints [39]. Dynamic postprocessing kinematic (PPK) technology can rapidly and conveniently link the camera's real-time positioning information with the ground GPS base station, thereby achieving coordinate correction and spatial interpolation for dense point clouds to obtain real geographic coordinate point clouds and DEMs. Currently, the highresolution topographic data obtained by SfM photogrammetry provides powerful basic data for geoscience research [40,41].
In recent years, the rapid development of low-altitude remote sensing platforms such as UAV has provided a new method of data acquisition for SfM photogrammetry. Because of its low cost, high data acquisition efficiency, and simple and flexible methods, this method has been applied to many fields of geoscience research [42,43]. As UAV flying heights are close to the ground, the accuracy of the obtained data is also improved, which can increase the terrain details and provide a new method for describing the tectonic active landscape [38,44]. SfM photogrammetry is also widely used to obtain topographic data of fault zones and facilitate the investigation of earthquake rupture zones and better describe earthquake disasters [45].

Acquisition of High-Resolution Topographic Data.
In this study, we used a fixed-wing UAV equipped with a fixed-5 Lithosphere focus digital camera to collect image data along the surface rupture of the 1931 Fuyun earthquake. The camera is a rigorously calibrated SONY Alpha7R full-frame digital measurement camera with 35 mm focal length lens and highquality 36 million-pixel full-picture image quality imaging, ensuring the data accuracy and color restoration of aerial survey images, which is conducive to the extraction and matching of later image features. The airborne platform is a fixed-wing UAV equipped with both a Global Navigation Satellite System (GNSS) and Inertial Measurement Unit (IMU) to ensure the stability of the flight route and airborne platform. Due to the integration of automatic cruise technology, the system can collect aerial photographic data according to the parameters set by the user (including flight route, flight altitude, speed, and image overlap).
The photographs were acquired on September 30, 2017, to October 11, 2017, in the late autumn, with fine weather and sparse vegetation on the ground, which was convenient for aerial photogrammetry (Figure 4). The flying height of the UAV was set to 300-350 m, with 85% and 70% forward and side overlap, respectively, and the image data were stored in JPG format. During the aerial image acquisition process, we used dynamic postprocessing kinematic (PPK) technology to determine the position of aerial imagery. Compared with the traditional image correction method based on ground control points (GCPs), this method can save a lot of time and labor costs. Furthermore, the plane positioning accuracy and elevation accuracy of the PPK method can reach centimeter level, which meets the accuracy requirements of data acquisition in this study. A total of five GPS base stations were installed along the aerial survey area, and the coordinates of the base stations were accurately located. The ground base station was placed in an area with a wide view and flat terrain, and the ground elevation angle of the satellite was greater than 15°. In the process of UAV aerial photography, the ground base station and the mobile station carried by the UAV observe the GPS satellite at the same time. During postprocessing, the GPS data received simultaneously from the base station and mobile station were linearly solved to form a virtual carrier phase observation for determining the relative position between stations. By setting the known coordinates of the base station, the aerial image coordinates were obtained. With the above parameters, the total number of collected photographs was 24,000, covering an area of approximately 110 km × 1:2 km.
After image acquisition, the aerial data were processed using the Agisoft PhotoScan software to register and construct the three-dimensional surface model. The data processing flow and parameter setting were performed according to the software documentation and other research papers [40,46]. Firstly, we inspected the aerial image visually to eliminate blurred and poorly exposed photographs. Secondly, the relative position of the aerial photographs was restored using the camera coordinates and image features for feature matching and tracking. Then, dense point cloud data with density > 50 m 2 were obtained by searching and matching between images according to the MVS principle. Finally, the Surfer software was used to interpolate the point cloud data by natural proximity interpolation to generate a DEM with the geographic coordinate system and a grid size of 1 m.

Identification and Measurement of Offset Markers.
When an earthquake occurs, the surface topography is distorted by fault activity to form fault landform markers, such as fault scarps and deflected gullies. These markers can be retained for a long time under suitable external conditions, thereby recording the displacements of multiple earthquake events. A basic assumption when using geomorphic markers to reveal the rupture history along a fault is that the occurrence rate of geomorphic markers used to measure displacement is (significantly) greater than the recurrence rate of (surface-rupturing) earthquakes [8]. Based on this 6 Lithosphere assumption, a suite of geomorphic assemblages formed between successive earthquakes will move together during subsequent earthquakes, recording seismic slip values. Tectonic and climatic changes can promote the formation of new geomorphic markers on a certain spatial scale [1,8]. On the one hand, a large number of well-preserved gullies are widely distributed in the southern foot of Altay Mountains. On the other hand, compared with the deformation rate of 10 mm/yr in the Tien Shan, the deformation rate of 2-3 mm/yr in the northeastern part of the Junggar Basin area and the long earthquake recurrence period of the Fuyun fault for nearly 10,000 years create favorable conditions for the formation of new geomorphic markers in interseismic periods [47,48]. Therefore, we believe that this region meets the above assumptions.
On a fault, the amount of slip formed by an earthquake will cause the offsets of geomorphic markers to generate clusters around some discrete values [49]. Therefore, multiple clusters represent the cumulative slip of a number of events, of which the cluster of the smallest displacement represents the displacement of the most recent earthquake [8,18]. Thus, we obtained the slip distribution of the Fuyun fault by measuring the horizontal distance of the staggered gully in the study area. The Dispcalc software based on Matlab was used to measure the horizontal displacement of the dislocation gully [11]. Firstly, this software determines the rupture position and the gully trace on the DEM, generating a longitudinal profile. Then, it uses "back-slip" to determine the horizontal displacement. Through multiple matching and fitting of the gully profile, the distance of horizontal movement is the horizontal offset of this gully when the cumulative probability is maximum. Finally, the software also provides the corresponding displacement back-slip map and confidence interval ( Figure 5).

Surface Rupture and Offset Observations. The 1931
Fuyun earthquake was characterized by a seismic surface rupture zone. Based on high-resolution topographic data and a field investigation, we mapped the surface trace and right-lateral offset of the gully. The linear characteristics of the surface rupture zone are obvious, and the geometric structure is almost continuous (Figure 6). The epicenter of this earthquake is located in Karaxingar. From the epicenter, the surface rupture zone extends to the north and south and 7 Lithosphere exhibits different features. The surface rupture zone to the north is straight, passing through Koktokay and extending to Kayleite. Conversely, the rupture zone exhibits many changes to the south, from NNW to nearly NNW, then to NNW, and finally to NNW. One branch turns toward a NW orientation and extends along the hillside, and another gradually disappears along the alluvial plain in front of the mountain. Moreover, some surface deformation S1 S2 S3 S4 S1 S2 S3 S4     Lithosphere characteristics are observed in this area, such as compression uplift and pull-apart basins, which can reflect the spatial geometric distribution of surface rupture. According to the geometric and tectonic geomorphic features, the surface fracture zone in the study area is divided into sections S1, S2, S3, and S4 from north to south ( Figure 6).
In this study, a total of 550 linear gullies with varying degrees of right-lateral dislocation were identified along the surface rupture zone, whose maximum, optimum, and minimum offsets were measured and recorded. The offsets, which vary from 0.9 m to 94.8 m, are plotted against their distance along the fault strike in Figure 7.

Offset Distribution and
Clustering. In order to analyze the large number of offset data, we calculated the Gaussian probability density function (PDF) for each offset measurement and differentiated offset clusters by summing each PDF to form a cumulative offset probability density (COPD) [50]. The values at each peak of the clusters represent the most frequently occurring values in the data set, and the width of the cluster reflects the discrete degree of the cluster. As most of the offsets (>90%) are less than 40 m, offset values of less than 40 m were used to calculate the COPD.
The rupture ranges of different earthquakes are not necessarily the same; therefore, it is more reasonable to analyze each rupture segment individually [50][51][52]. The offset distribution of the four rupture segments falls into multiple separate clusters (Figure 8). S1 offsets ranged from 2.

Data Comparison.
The collection of satellite remote sensing image data is easily affected by weather factors such as snow and clouds, and the corresponding elevation information is lacking if it is not corrected. Conversely, the dense point cloud data obtained by the SfM method can generate three-dimensional and terrain data, which color attribute makes the point cloud data more stereoscopic and realistic; therefore, the obtained texture features can be more easily applied to visual interpretation [53], which is beneficial for detailed local geomorphology research. In order to evaluate the reliability of generating terrain data from SfM photogrammetry, we compared the ability of the DEM and Quickbird satellite images to interpret gully morphology in the same location [20]. In Figure 9, the yellow lines are the interpreted gullies. The results show that the terrain data generated by the SfM method have less uncertainty, which enables us to fully capture the morphological characteristics of gullies and identify previously undiscovered gullies (Figure 9(d)), thereby improving the accuracy of horizontal offset measurement.
We also compare the offset results of 54 groups with the offset data measured by Klinger et al. [20] in the same location ( Figure 10). The R 2 value of the fitting degree between the two data types is 0.002. Most of the offset values differ 10 Lithosphere by 1-2 m; however, some of the offset values can be as much as three times the measured values. In addition to subjective factors for determining the shape of the gully, previous offset data may record the cumulative value of the horizontal offset from multiple earthquake events. Figure 11 is a good example. This place is a piedmont alluvial fan, on which a series of fault gullies provide us with a lot of new offset data. The gully in the white frame is measured, and the result is 6.7 m (blue line). And in the downstream of the gully, we find a gully with undercutting marks but not obvious along the sliding direction of the fault (yellow line). We think it should have been the downstream of this gully, but it has been abandoned and strongly modified due to the movement of the fault. The measured distance from the upstream gully is 13.7 m, so we believe that this result should be a cumulative offset. This is also one of the reasons for the big difference between the two groups of data.
In order to evaluate the accuracy of this offset measurement, we compared our offset measurement results with field measurements at the same location. Forty dislocation trenches with good linear traces are selected for field mea-surement, and the results show good correlation between the two data ( Figure 12), with a linear fitting equation of f ðxÞ = 0:805x + 1:21 and an R 2 value of 0.837. The difference between the two data sets may be caused by the shape of the dislocation trench, the fitting method, and the measurement error itself. Error in the measurement is mainly caused by uncertainty of the fault strike, and determining the gully shape as morphological analysis is influenced by subjective factors. In addition, if different fault strikes are selected, the distance between the trenches will also change. The part of the gully near the fault is affected by the fault activity multiple times, which greatly affects its shape; therefore, it is not suitable to measure the distance in this location. Therefore, when we use the Dispcalc software, one person was chosen to match and fit the fault and gully trace 6-10 times for each measurement, selecting the relatively straight and stable part of the gully trace to reduce uncertainty in the measurement.  Figure 11: Different offsets are obtained by measuring the same gully. Normal value is the offset obtained from this study, and blue italic value is the offset measured by previous study [20]. 11 Lithosphere distribution of seismic events, thereby revealing the rupture history and recurrence patterns of strong earthquakes and providing valuable information about the potential future behavior of faults [54,55]. Lin et al. [52] concluded that the interpretation of paleoearthquake events by COPD should ensure the following three points: (1) original historical coseismic slip data are used; (2) the coefficient of variation (CoV) value is relatively low; (3) there are sufficient offset measurement data. In this study, the coseismic offset of the Fuyun fault caused by the 1931 earthquake is well recorded, the CoV value is used to evaluate the regularity of cumulative offsets, and the global historical coseismic surface rupture segments have an average CoV value of approximately 0.58 [52,56]. Here, we follow the method described in previous studies to calculate CoV values of 0.23, 0.26, 0.32, and 0.28 for each rupture segment [8,18]. Moreover, the amount of measured data is sufficient to meet the above requirements. Therefore, we believe that the paleoearthquake events of the Fuyun fault can be explained by the COPD.
In the offset distribution plot, the clustering characteristics of segments S1 to S4 show that the average offsets of the most recent events in each segment are 6.0 m, 5.3 m, 5.1 m, and 6.1 m, respectively. Within section S1, the cumulative offsets of 9. mately two and four times as large as the coseismic offset of 6.1 m for the most recent earthquake event. Thus, there is a clear relationship between the cumulative offsets and the most recent event in each segment of the fault, whereby the cumulative offsets of two, three, and four times may correspond to the cumulative offsets of two, three, and four earthquake events, respectively. This indicates that, in addition to the 1931 earthquake, there have been three historical earthquake events of the same magnitude. Although no cluster corresponding to three times the coseismic offset is found in segment S4, the cluster width corresponding to two and four times the coseismic offset is obviously larger than that of the other three segments, indicating greater discretization of the two clusters, that is, a larger range of data comprising the cluster. This may be due to the fact that, in the second to most recent event, the small value of the offset belongs to the most recent cluster, and the large value belongs to the fourth from last cluster, which also explains why the cluster corresponding to three times the coseismic offset does not exist. According to the offset distribution of each segment, the coseismic offsets larger than 25 m do not form obvious clusters; therefore, these results identify four earthquake events at most but are unable to identify five or six events like previous studies [20].
The clusters of minimum offsets in adjacent segments are connected to establish the coseismic offset distribution of the 1931 earthquake ( Figure 13). The results show that the coseismic offset of the 1931 earthquake varies from high to low along the fault strike; the coseismic offset increases rapidly at 0-5 km, 28-33 km, 48-52 km, and 74-76 km and decreases rapidly at 26-28 km, 47-48 km, 70-71 km, and 105-110 km. In this way, four connected waves are formed, with peaks at 18 km, 36 km, 53 km, and 71 km, respectively. There is a high gradient of offset on both sides of the wave, and three grooves are formed between adjacent waves at

Rupture Model of the 1931 Fuyun Earthquake.
Understanding the rupture model is a necessary condition for earthquake risk assessment. In an earthquake, the process and propagation of rupture are usually limited by the geometric complexities of the fault, such as bends or step-over areas. These complex zones divide the fault into segments [57]. A single segment can be ruptured independently in different earthquakes, or multiple segments can participate to create a cascade rupture [58]. Whether the earthquake rupture can pass through the discontinuous zone is controlled by two main factors: the offset and offset gradient when the rupture enters the step-over zone and the width or bend angle of the step-over area. A rupture with a higher offset gradient has a greater possibility of breaking through the discontinuous zone, whereas a lower offset gradient has less potential to break through the discontinuous zone [59]. At the end of the rupture segment, the coseismic offset distribution usually takes a "rainbow" form ( Figure 14), reflecting the stress drop and characteristic offset distribution of a complete fault segment [60]. As for the width of the stepover area, a combined analysis of historical earthquakes and numerical simulations showed that the critical width of a rupture terminating in the step-over area is 3-4 km; thus, historical earthquake ruptures wider than 3-4 km cannot continue to expand through the discontinuous step region [61]. Moreover, Lozos et al. [62] simulated the critical bend angle of a compressional bend rupture and showed that a compressional bend angle of less than 18°can make the rupture propagate through the entire fault. Therefore, we discuss the possible rupture patterns of earthquakes  Figure 14: Schematic of the two rupture termination modes: "dogtail" and "rainbow" (modified from Ward [60]). The fracture extending to the right encounters a fault segment with the strength of 90 bars and an initial stress of 25 bars. This "hard" barrier prevents the rupture from spreading within 1-2 km of the boundary of the segment, forming a "rainbow" sliding pattern. The expansion fracture on the left side occurs on a segment with uniform strength of 25 bars, but only encounters preexisting stress of 5 bars, forming a "dogtail" like sliding pattern.  13 Lithosphere according to the coseismic offset distribution and surface rupture geometry of the 1931 Fuyun earthquake.
In the offset distribution map of the 1931 Fuyun earthquake, the offset shape of each segment is convex and has a high offset gradient at the boundary of the step area ( Figure 13). However, there is one exception: at the northern end of section S1, the offset only decreases from 6 m to 4 m, indicating that the boundary of this section is larger than this value. The field investigation confirmed this finding, indicating sporadic surface ruptures in the northward area of the surface rupture ( Figure 15) where there should be a smaller coseismic offset, resulting in the rapid decline in the offset value of the S1 section boundary. This form of offset distribution reflects four complete stress drops, indicating that the surface rupture in each segment represents a completely independent fault, and that surface rupture is likely to break through the step-over propagation. Mapping of the rupture of the 1931 Fuyun earthquake confirms that the seismic surface rupture zone is composed of segments separated by several step-overs. Specifically, there is 1.5 km wide compressional step-over zone between S1 and S2 and a 0.8 km wide step-over zone between S3 and S4. These step-overs are not wide enough for the rupture to end here. Regarding the bend at the junction of S2 and S3, the measured bending angle is 17°. Bending of this scale does not exceed the results simulated by Lozos et al. [62], indicating that the rupture can also propagate here. Based on the above analysis, we suggest that the surface rupture of the 1931 Fuyun earthquake was a cascade rupture formed by four rupture segments.

Conclusion
In this study, we provide new mapping results of the surface rupture of the 1931 Fuyun earthquake using UAV-based SfM photogrammetry. Through detailed mapping of the topographic features of the surface rupture zone and measurement of 550 groups of fault gullies, we obtained the surface morphology of the Fuyun fault in higher resolution than previous studies, revealing four distinct rupture segments. We also obtained a more complete distribution of slip offsets and identified four earthquake events from the clustering features, including the 1931 earthquake. Therefore, this study does not support the finding of previous studies, which identified five to six earthquake events [20]. By extracting the offset profile of the 1931 earthquake, we revealed that the slip distribution along the fault strike shows four continuous undulations, corresponding to the four surface rupture segments, and a high offset gradient in the step-over area connecting the rupture segments. Combined with the width and bend angle of the step-over area at the rupture segment boundaries, this indicates that the 1931 Fuyun earthquake was a cascade rupture formed by four distinct rupture segments.
This research confirms the application value of highresolution topographic data obtained by UAV-based SfM photogrammetry for detailed geological mapping and fault slip distribution analysis. This technique enables us to collect topographic and geomorphic data more quickly after an earthquake, provide basic data support for postearthquake disaster relief work, and greatly reduce manpower and material costs. Furthermore, this method can more effectively restore the slip distribution and seismic propagation characteristics of faults, thereby providing valuable data for seismic risk assessments.

Data Availability
All gully displacement data can be found in the supplementary file.