Similarity Analysis: Revealing the Regional Difference in Geomorphic Development in Areas with High and Coarse Sediment Yield of the Loess Plateau in China

: The development of loess landforms is controlled by underlying, pre-existing paleotopography. Previous studies have focused on the inheritance of loess landform and the control of underlying paleotopography on modern terrain based on the digital elevation model (DEM), while the similarities and differences between modern terrain and underlying paleotophotography have not been directly spatialized. In this study, areas with high and coarse sediment yield (AHCSY) in the Loess Plateau of China were selected to form the study area, and the DEM of the study area’s underlying paleotophotography was reconstructed using detailed geological maps, loess thickness maps, and underlying paleotopographic information. The hypsometric integral ( HI ) and spatial similarity analysis methods were used to compare the spatialized difference between underlying and modern terrain of the Loess Plateau from the perspectives of the landform development stage and surface elevation, respectively. The results of the HI method demonstrate that essentially, there are similarities between the geomorphologic development stages of underlying and modern terrain, and only some local differences exist in some special areas. The results regarding the spatialized coefﬁcient of relative difference and the Jensen–Shannon divergence demonstrate that the thicker the loess is, the weaker the similarity is, and vice versa. Meanwhile, according to the present loess landform division, the order of regional similarity from low to high is as follows: loess tableland, broken loess tableland, hilly regions, dunes, and the Yellow River Trunk. The use of the similarity analysis method to analyze similarities between underlying and modern terrain plays an important role in revealing the inheritance of loess landforms.


Introduction
The Loess Plateau of China is the region with the most concentrated loess distribution and the deepest covering thickness in the world, and it is also one of the regions with the most serious soil erosion [1,2]. The formation of the Loess Plateau has experienced a long process of erosion, denudation, and accumulation. A large amount of dust accumulation has covered the underlying terrain since before the Quaternary, forming a large area of dust accumulation [2]. Although the deposition of loess has a certain influence on the shaping of the topography of the region, it can only be said to be a minor modification of the original underlying terrain [3]. Therefore, there is a certain similarity in space between the modern and underlying terrain before loess deposition. The study of the spatial similarity and differences in distribution characteristics between the paleotopography before loess deposition and the modern terrain is significant for understanding the trend and situations that surround loess development [4].
The Loess Plateau, as a key area of geomorphological research, has always attracted attention. A large number of scholars have conducted qualitative and quantitative research on the erosion state and development stage of the loess landform [5][6][7]. A great amount of the terrain of the modern Loess Plateau has inherited the complex and diverse formats that were present before the Quaternary [2]. Most loess hills on the Loess Plateau are controlled by the underlying paleotopography; primary loess hills are mainly formed by loess accretion on the basis of the underlying paleotopography, and secondary loess hills are transformed by gully erosion [8]. Li and Lu [9] estimated the average surface erosion rates on the Loess Plateau in each glacial or interglacial period over the past 250 ka as well, and the primary results demonstrate that relatively intense erosion has existed on the Loess Plateau in the recent geological past. Based on research results regarding the variation characteristics of the strata structure and paleotopography in the north and south of the loess in the northern Shaanxi, Guo [10] proposed and demonstrated that there was an unconformability between the loess landform (tableland, beam, and hill) and underlying paleotopography in northern Shaanxi. A systematic particle size analysis of the Miocene loess profile demonstrated that Miocene and Quaternary loess and paleosol samples had similar characteristics in particle size distribution in Qin'an, Gansu [11]. From this research, it can be observed that the relationship between modern and underlying terrain needs to be further researched. In recent years, with the rapid advancement of data acquisition methods and the improvements in digital terrain analysis methods, many achievements have been made in the study of loess landform through digital terrain analysis [12,13]. Chen et al. [14] established a digital model of ancient landforms in the Jurassic and Tertiary periods in a small river basin in the Fugu county of Shaanxi Province, China, based on the application of high-precision GPS, RS, and GIS. It was found that the modern terrain in the river basin inherited the erosion surface present in the Jurassic and Tertiary periods to a large extent. Zhu et al. [15] studied the HI and its spatial variation in the Loess Plateau. Based on the DEM, Cai et al. [16] studied the fractal characteristics of typical geomorphic forms on the Loess Plateau in northern Shaanxi. Research regarding the inheritance on the landform evolution of the key soil erosion areas of the Loess Plateau demonstrated that the underlying paleotopography of the pre-Quaternary bedrock controls the evolution of the loess-covered landform to a large extent [17,18].
Previous studies have mainly focused on the controlling effect of the original, underlying pre-Quaternary surface layer, before the evolution of loess deposition on the present landform [18], and the inheritance analysis of the loess landform [19]. However, little attention has been paid to the spatial similarity between the modern terrain and underlying paleotopography. Similarity analysis is one of the most widely used methods in geographical research [20]. The morphology and distribution of the paleotopography may have significantly affected spatial variation in the loess deposition process [18]. The analysis of the similarity between the modern and underlying terrain on the Loess Plateau is helpful for understanding the development trend of the loess landform, especially in determining whether there are great differences in the similarity between two landforms in different loess landform areas. At the same time, it is helpful to more intuitively understand the spatial distribution difference between the underlying terrain and the modern terrain, which can also intuitively highlight the strong control effect of the underlying terrain and the weak control effect of the underlying terrain from the space, so as to help further analyze which regions may be inherited gullies and which regions are modern erosion gullies. This paper selected areas with high and coarse sediment yield (AHCSY) in the middle reaches of the Yellow River as the study area. Based on the geological map data, the early measured loess thickness distribution map, high-resolution remote sensing images and DEM data, we constructed a DEM of the underlying paleogeomorphology in the study area. The spatialized coefficient of relative difference and the Jensen-Shannon (JS) divergence reveal the spatial similarity distribution characteristics of the Loess Plateau and paleogeomorphology on a macro scale.

Study Area
The study area (Figure 1) was located in areas with high and coarse sediment yield (AHCSY) in the middle reaches of the Yellow River, which is one of the areas with the most serious soil erosion and the most fragile ecological environment in China [21,22]. It is located between Hekou Town and Longmen District of the Yellow River and the upper reaches of Jinghe and Luohe River, and the main landform types are loess hilly gullies and loess tableland gullies [23]. Meanwhile, the gully density is large and the slope is steep and deep, and the surface soil in the area is mainly dominated by erosive sandy loess. It is in an arid and semi-arid area with little and very concentrated rainfall [24,25]. The average annual precipitation is about 450 mm, decreasing from southeast to northwest. More than 60% of the annual precipitation occurs from July to September and is often in the form of heavy rain [25]. By 2020, the area of soil erosion in the middle reaches of the Yellow River reached 107,980 km 2 , accounting for 54.33% of the total area. The intensity of soil erosion is mainly mild and moderate erosion, of which 54% is mild erosion, 28% is moderate erosion, and 18% is strong and above erosion [25,26].

Material and Underlying Terrain Modeling
In this study, we used 90 m resolution data from the Shuttle Radar Topography Mission (SRTM) as the basic elevation information source, and its accuracy was sufficient to ensure the expression and mapping of large area scales [4,27]. The geological map of the Yellow River Basin and the loess thickness distribution map of the Yellow River Basin served as the basic data regarding loess thickness [21]. Firstly, using the geological map of the Yellow River Basin (Figure 2), the bedrock outcropping points of different geological ages and rock types marked in the study area were found, their coordinate positions and altitudes were collected, and the high-definition remote sensing imagery map (the 2.5 m resolution data of Map World) was applied for georeferencing to locate the obtained bedrock outcropping points. In this way, we obtained all of the bedrock outcropping points in the study area. Secondly, the contour line of the loess thickness was obtained from the loess thickness map of the Yellow River Basin. The distribution points of the loess thickness were obtained by converting the contour lines of the loess thickness into points. Then, the actual measurement data and loess thickness contour map created by Liu [3] were used to encrypt the loess thickness points. Finally, we obtained the sampling points of loess thickness distribution in the study area, shown in Figure 3. Based on those sampling point data, we constructed a digital elevation model of loess thickness using interpolation. In the selection of the interpolation method [28,29], we chose the spline function (root mean square error (RMSE) 19.8 m) because of its high accuracy and relatively low variation [30]. Using a simple mathematical operation for modern DEM and the loess thickness model, the height of the modern DEM minus that of the loess thickness, the paleotopographic surface model could be determined for the study area ( Figure 4).

Hypsometric Integral
Hypsometric integral (HI) is an important indicator that reveals watershed geomorphic morphology and development characteristics through statistics of watershed surface elevation combination information [5,6]. There are many ways to calculate the HI value [31]. Here, we mainly adopted two methods to calculate the HI value: the hypsometric curve method and the elevation-relief method. Using the hypsometric curve method, the HI value is calculated by its definition [32], which refers to the relationship curve between the horizontal section area of the experimental sample area and its relative elevation above the estuary. The rectangle, whose vertices are the upper and lower endpoints of the hypsometric curve, is divided into two parts by the elevation curve. The ratio of the area under the curve to the total area of the rectangle is the HI value of the watershed. The formula is as follows: where HI is the value of hypsometric integral. H is the watershed topography difference, A is the watershed area, h is the relative elevation of the contour line, and a is the horizontal section area cut by the contour line. The elevation-relief method is the most efficient and simple method in calculating the HI value. It is a method of estimating the HI value derived by Pike and Wilson through mathematical formulas [33]. The calculation formula is as follows: where HI is the value of the hypsometric integral. H mean indicates the average elevation, H min indicates the minimum elevation, and H max indicates the maximum elevation. In the paper, formula 2 was used to calculate the HI of paleotopography and the modern terrain. The selected tool was Zonal Statistic in ArcGIS.

Spatialized Coefficient of Relative Difference
In previous studies, the coefficient of relative difference is a quantitative method that can only provide a value of statistical similarity between two layers, such as 0 or 0.5. In this study, we improve this method so that it can express local spatial similarity at different neighborhood scales, such as a neighborhood of 15 × 15, so we call it the spatialized coefficient of relative difference [34]. The spatialized coefficient of relative difference is used to quantitatively analyze the similarity between the modern terrain and paleotopographic DEM. The DEM layers of the modern terrain and paleotopography are A i base and A i comp , respectively, and the relative difference coefficient α can be defined as [34]: where A i base represents the elevation value of the modern terrain at the i-th grid unit, A i base is the mean value of elevation in a certain spatial extent, and A i comp is the elevation value of the underlying terrain at the corresponding position. The relative difference coefficient, α, describes the overall deviation between the modern terrain and the underlying terrain. If α = 1, the two DEM layers are exactly the same. The smaller the α is, the greater the difference between the two DEM layers. When α < 0, the two DEM layers are the same.

Jensen-Shannon Divergence
The Jensen-Shannon (JS) divergence, as an algorithm for measuring similarity, is often used in landscape similarity searches [35,36] and aerial remote sensing image change monitoring [37,38]. Here, we used use the JS divergence to analyze the similarity between modern terrain and underlying paleotopographic DEM. The JS divergence is derived from the Kullback-Leibler (KL) divergence. The definition of KL divergence is as follows [39]: As KL divergence is asymmetric, that is, KL( f (x) g(x)) = KL(g(x) f (x)), in order to make it symmetric, Jansen and Shannon proposed a new method for calculating the relative entropy, that is, taking the average of both sides of the inequality [40,41]. Namely: In the formula, f (x) represents the DEM of the paleotopographic surface, and g(x) represents the DEM of the modern surface. The value of JS divergence is [0, 1]. When the value of JS is 0, the two images are completely the same. The closer the JS value is to 1, the worse is the similarity. Figure 5 shows the difference in the hypsometric curve between the modern and underlying terrain. The HI value of the modern terrain is about 0.51, while the HI value of the underlying terrain is about 0.49. The hypsometric curve of the land surface parallelly rises, and the trend of hypsometric curve is very similar, which conforms to the inheritance characteristics in the process of loess deposition [4]. The spatial distribution diagram of HI values of the modern and underlying terrain is shown in Figure 6. At the same time, according to the spatial distribution map, we calculated the area proportion of the modern and underlying terrain in different development stages of the whole study area. As shown in Table 1, which displays the area proportion of different geomorphologic development stages, there are some slight changes in the area ratio between the stages of geomorphologic development of modern terrain and that of underlying terrain. As can be observed from the spatial distribution map of the HI values of the modern and underlying terrain, in most areas, the HI value is between 0.35 and 0.6, in the mature stage of geomorphologic development. The HI value of the Yellow River Trunk and the outcropping of the bedrock is less than 0.35, and this area is in the last stage of geomorphologic development. In areas of intense wind erosion near Ordos city, the HI value is greater than 0.6, and the landform development is in its juvenile stage. In the loess tableland and broken loess tableland centered around the Huanxian county, and the loess hilly area centered on the Wu Qi county, the HI value of the underlying terrain is mostly 0.35~0.43, while the value of modern terrain is higher, at 0.43~0.52, indicating that the modern terrain is more eroded than the underlying terrain. According to the research results of Zhao et al. [42], the stage of geomorphologic development has a certain relationship with the erosion intensity. In areas with severe erosion, most of the landforms are in the juvenile stage of geomorphologic development, the areas with moderate erosion are in the mature stage, and the areas with mild erosion are in the late stage [43]. At the same time, near the Yellow River Trunk and other bedrock channels, most of the modern terrain is in the mature stage of geomorphologic development, indicating that modern erosion is intense and hydraulic erosion is still increasing [42,44,45].

Similarity Analysis
The coefficient of relative difference and JS divergence are both used to measure the similarity. In previous studies, they were presented in the form of numerical values [34,35,37,41]. In this paper, the coefficient of relative difference and JS divergence were spatialized to measure the spatial similarity of the two DEM layers. The closer the value of the coefficient of relative difference is to 1, the higher the similarity, and the lower the similarity is when it is closer to 0; the closer the JS divergence is to 0, the higher the similarity is, and vice versa, the lower is the similarity.

Spatialized Coefficient of Relative Difference
When we calculated the coefficient of the relative difference of the two maps in a statistical way, the value of α was 0.84. Through this value, we could only determine that the underlying and modern terrain have high similarity; we could not know their spatial similarity distribution. However, through the calculation of the spatialized coefficient of relative difference, we could intuitively see the spatial distribution differences between the two layers. When using the spatialized coefficient of the relative difference to calculate, the cell size of the neighborhood analysis will have a certain impact on the result. We compared and analyzed the calculation results of different neighborhood analysis cell size (Figure 7) and made statistics on their characteristic values ( Table 2). The cell size was between 27 × 27 and 42 × 42, the calculated coefficient of relative difference was stable, and the calculated values beyond this range have certain errors. Therefore, we compromised and selected the cell size of 33 × 33, as shown in Figure 8. As can be observed from the figure, in the northern part of the study area, the Yellow River Trunk and the bedrock outcrop area, such as Qingshuihe county, Hequ county, Baode county and Wuqi county, the relative difference coefficient value is close to 1. This demonstrates that the underlying terrain in these areas have high spatial similarity with the modern terrain. However, in the central and western regions of the study area, the values of the coefficient of the relative difference are close to 0, indicating that the modern terrain and the underlying terrain are not similar.

Jensen-Shannon Divergence
The spatial distribution of the JS divergence on the modern terrain and the underlying terrain is shown in Figure 9. By comparing the two figures in Figures 8 and 9, it can be observed that the spatialized coefficient of the relative difference and JS divergence are consistent on the whole, but there are some local differences. For example, in the windy and sandy area near Ordos, the spatial similarity calculated using the two maps is different, the relative difference coefficient results are not similar, and the JS divergence results are similar. Perhaps because neighborhood analysis is not involved in the calculation of JS divergence, the result is relatively smooth. In the western region with the Huanxian county as the center, the JS divergence is close to 1, and the similarity is low. In the central and eastern regions, the value of JS divergence is between 0 and 1, showing partial similarity, for example, in the Shenmu area, Suide county, Lvliang area, Yan'an area, and Huachi county. At the same time, through the comparison with the loess thickness map (Figure 2), we know that the distribution of similarity is related to the thickness of loess accumulation to some extent-the thicker the loess deposits, the worse the spatial similarity between the modern and the underlying terrain. For example, the thickness of loess in the surrounding area of the Huanxian county is close to 300 m, the coefficient of relative difference is close to 0, the deviation of JS divergence is the largest, and the similarity is weak. This indicates that the influence of underlying terrain is stronger in the region with thin loess coverage than in the region with thick loess accumulation. In the early stage of loess accumulation, the underlying terrain plays a decisive role, but in the area of thick loess accumulation, its influence is relatively weak [19]. The loess in China is generally believed to be transported and accumulated by wind, and the thickness of the loess tends to decrease from west to east and from north to south [1]. The study area was located in the Cenozoic tectonic movement area of the Ordos stable block. At the beginning of the Quaternary, as the global climate became colder and the East Asian monsoon strengthened, the Ordos block received more dust accumulation. Moreover, the Quaternary desert appeared in the northwest, and a large amount of dust accumulation covered the underlying terrain before the Quaternary [46]. During the development and evolution of the loess, it was damaged by various types of natural erosion (hydraulic erosion, gravity erosion, and wind erosion) and human activities [47][48][49][50], resulting in the fragmentation of the surface and the development of gullies. Tableland, beam, and hill are the three most common landforms in the Loess Plateau [51]. There are also great differences in similarity between paleotopography and modern topography under different geomorphic forms. We overlaid the geomorphological map, and the spatial distribution map of similarity are overlaid; the result is shown in Figure 10. At the same time, we calculated the JS divergence of typical landforms in the study area, and used the boxplot to show the variation of the JS divergence of different landforms, as shown in Figure 11. It can be observed from the figure that the underlying terrain of the loess tableland and the broken loess tableland has a weak similarity with the modern surface, which is specifically manifested in the Huanxian county and other areas. In the loess tableland and the broken loess tableland, the similarity between paleotopography and modern terrains is weak, for example, in Zichang, Yan'an, Suide, and other areas. In the dune area and the Yellow River Trunk, the similarity between the underlying terrains and the modern terrains is relatively high, for example, in Yanchuan, Baode, Hequ, Qingshuihe and other areas. It can be observed that the relief of the modern terrain is roughly similar to that of the paleosol layer. The paleosol layer is relatively flat in the tableland area, while it is inclined to the direction of the adjacent large gullies in the ridge area [52].

Terrain Profiles Analysis
In order to more directly reflect the spatial relationship and distribution pattern between the paleotopographic DEM and the modern surface DEM, we obtained a set of nests of the profile line of the paleotopographic DEM and the modern surface DEM in the study area. As shown in Figure 12a, eight pairs of terrain profiles were extracted in the whole study area, lines a, c, and d from east to west, lines b and e from south to north, lines g and h from northwest to southeast, and line h from northeast to southwest. Figure 13 shows the altitude graphs of the terrain profiles extracted from the paleotopographic DEM and the modern surface DEM. Combined with the geomorphic map, profiles a and b in Figure 12a are mainly located in the loess tableland and loess broken tableland; the other terrain profiles are mainly located in the loess hilly and gully regions. Profiles d, f, and h crossed the Yellow River, and profiles c, e and g crossed the tributaries of the Yellow River, including the Wuding River and Yanhe River. At the same time, we extracted the JS divergence of each profile line and drew a boxplot, as shown in Figure 12b. As can be observed in the eight pairs of terrain profiles, there are great differences between ancient and modern terrains in the loess tableland and loess broken tableland, while the similarity is high in hilly and gully regions. Meanwhile, in the macroscale of the study area, a significant landform inheritance relationship can be observed between these two terrains during the loess deposition process [19].

Evaluation of the Loess Thickness
In this paper, the loess thickness is simulated by the interpolation method, and the surface of the paleotopography in areas with high and coarse sediment yield of the Loess Plateau is reconstructed. In the selection of the interpolation method, we referred to the interpolation method adopted by Xiong et al. in the reconstruction of the underlying terrain, and selected the spline function with high precision and small variation [19], with an RMSE of 19.8 m. The error range was reasonable in the study of using DEM to reconstruct the paleotopography [4,18,19]. For example, Xiong et al. used the spline function interpolation to reconstruct the paleotopography, and the RMSE was 36.5 m [4]. In addition, there are other methods for simulating loess thickness. For example, Zhu et al. generated the loess thickness map by using the digital elevation map (DEM) of the Loess Plateau and the neighborhood analysis in ArcGIS software [52]. In contrast to the interpolation method, this method is complex and depends on the size of the grid. Furthermore, due to the limited means of data acquisition, there is a difference between the calculated and the actual loess thickness. Through field investigation, it was found that in some areas the erosion base level is much lower than the actual base level, and even the rock is cut down by tens of meters, especially in the main gully. The soil thickness value may be much higher than the actual thickness value [53]. Near the Yanhe River Basin, Mizhi county, and Suide county, either the main river channel or a small branch ditch was cut below the bedrock, and the large branch ditch was even cut tens of meters deep. Therefore, it is necessary to obtain more accurate loess thickness data through extensive field profile survey and borehole sampling in the future.

Evaluation of the Similarity Analysis
Based on the improved similarity method, this study discussed the spatial similarity between the modern terrain and the paleotopography. Although the spatialized coefficient of relative difference and JS divergence are consistent on the whole, there are some local differences. When using the spatialized coefficient of relative difference to calculate, the cell size of the neighborhood analysis will have a certain impact on the result. However, in the process of calculating JS divergence, the influence of the analysis window is not involved, and the results obtained are relatively smooth. Therefore, when we use these two methods to calculate the spatial similarity of the two layers, JS divergence is a good choice. For the spatialized coefficient of relative difference, we should carefully consider the uncertainty of it and pay attention to the scale range of the analysis.
The similarity analysis results demonstrate that the thicker the loess is, the weaker the similarity is, which proves that the similarity between the modern terrain and the paleotopography is related to the thickness of the loess. This finding is consistent with that of Xiong et al. [4], who reported that the thicker the loess is, the less the underlying paleotopography controls the modern terrain. What is more, in the previous studies, through the comparative analysis of the paleotopographic morphology and the development history of individual loess profiles before the deposition and the loess profiles overlying them, it is found that the inheritance of the loess has regional characteristics. For example, Liu and Qiao et al. demonstrated that the loess landform had inheritance [2,11], while Guo's research demonstrated that there was no significant inheritance in the Luochuan Tableland and hilly regions [10]. This is consistent with the results of our similarity analysis, which demonstrate that the order of regional similarity from low to high is the loess tableland, broken loess tableland, hilly regions, dunes, and the Yellow River Trunk. In addition, the results of this study deepen the understanding of the development of gullies on the Loess Plateau, implying that most of the gullies in the loess tableland are modern erosion gullies, and most of the gullies in the hilly region are inherited gullies [54]. At the same time, although the underlying terrain has a profound influence on the formation of the modern loess landform in the overall pattern, the underlying terrain and modern erosion process have a joint influence on the development trend of the modern loess landform in the local landform object [42]. Therefore, the inheritance of gully erosion to paleo-eroded zones should be considered as an important factor in the study and prediction of gully erosion trends in the modern Loess Plateau [55]. Notably, although this study revealed the relationship between the modern terrain and the paleotopography at a macroscale, the results need to be further studied through geophysical ground data and borehole data.

Conclusions
The similarity analysis method used to study the similarities between the underlying and modern terrain plays an important role in revealing the inheritance of loess landforms. The proposed analysis methods, such as HI and the spatialized coefficient of relative difference and JS divergence, are effective in spatially displaying the changes in underlying and modern terrain in loess landforms, and the method is expected to be applied to the similarity analysis of loess landforms at other scales.
The results of the spatial similarity analysis demonstrated that there is a significant similarity between the modern and underlying terrain, which proves the inheritance relationship of loess landforms. The hypsometric curve for the modern terrain was demonstrated to be higher than and nearly parallel to the underlying terrain with a similar shape and trend; it was in line with the inheritance characteristics of the loess deposition process. When the HI value maps between paleotopography and modern terrain were compared, they were also shown to be essentially similar, with local difference only occurring in some special zones. The results of the similarity analysis demonstrated that the high spatial similarity between the modern terrain and paleotopography were concentrated in the north of the Yellow River Trunk. According to the division of loess landforms, the order of regional similarity from low to high is the loess tableland, broken loess tableland, hilly regions, dunes, and the Yellow River Trunk. The similarity distribution between the underlying paleotopography and the modern terrain were demonstrated to have a certain relationship with the distribution of loess thickness. The thicker the loess deposit accumulation is, the weaker is the similarity. The results also imply that most gullies in the loess tableland and the broken loess tableland may be modern erosion gullies, and most gullies in hilly areas may be inherited gullies. Although the Loess Plateau has been eroded to varying degrees in the process of development and evolution, its overall shape is similar to that of the underlying terrain.
The employed methods were found to be useful but could only be applied to macroscale analyses due to the limited data sampling density. The development of modern geophysical techniques would provide better conditions for the future study of loess landform evolution, because it will enable more precise investigation and modeling. In the future, a better sample area should be selected to build a digital elevation model of the loess surface, so as to deepen the understanding of the genesis and development mechanism of loess geomorphology.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.