Past rainfall-driven erosion on the Chinese loess plateau inferred from archaeological evidence from Wucheng City, Shanxi

The history of changes in the Earth’s surface can be investigated by numerical modelling of landform evolution. However, such models rely on evidence which is often removed through erosion or other surface processes. Here, we use archaeological observations from the remains of Wucheng City in Shanxi, China, which was built on a loess area during the Neolithic period, to reconstruct the paleo-surface of the area and inform a landform evolution model. We identify differences between the present-day and paleo-land surface which suggest the topography was previously much less rugged. We find that period of heightened erosion rate in the early Holocene coincides with the period of increased rainfall, suggesting climate was the primary control on landform evolution during this time. We conclude that paleo-surface observations taken from archaeological sites, as demonstrated by this study, are valuable tools to inform landform evolution models in other regions where evidence is limited. The palaeo-topography around the ruins of ancient Wucheng City, China suggests erosion rate during the Holocene was largely driven by changes in precipitation, and shows how archaeological evidence can inform landform evolution models.

T he Earth's surface consists of different landforms. These landforms act as dominant factors constraining the Earth's surface processes, which profoundly influence the climate mode, hydrological process, soil type, landscape ecology, biological diversity and human activities in an area [1][2][3][4][5] . As landforms change via Earth's interior and exterior forces, the related phenomena and processes coevolve during this process 6,7 . Thus, modelling the evolution of landforms is one of the most important aspects of understanding the causes and consequences of topographic change, as well as the genesis of the Earth 8,9 .
For a long time, the landform evolution process was understood via the geomorphic cycle model 10 . As geomorphology developed, the mechanism of landform evolution was also revealed through field work by investigating the material composition 11,12 , formation age 13,14 , and environmental background of an area 15,16 . These studies successfully helped to deepen the understanding of landform evolution from the perspective of formation mechanisms. In recent years, advanced earth observation technology has been used to quantitatively express surface morphologies from the macroscale to the microscale [17][18][19][20] . The landform evolution process has been further explored by using space-for-time substitution 21,22 or dynamic monitoring in a small area 23 . Currently, different numerical landform evolution models have been proposed to deepen the understanding of physical geomorphic processes and mechanisms 3,8,24 ; these models have helped in analysing how landform morphology responds to Earth's interior and exterior forces [25][26][27] . However, landform morphologies simulated by numerical models are mathematical simulations of surface morphologies; whether the simulated surface can be used as an actual paleo-surface still needs to be verified with evidence 28,29 . Thus, it is necessary to reconstruct a reliable paleo-surface to verify a simulated surface morphology, but evidence is rarely available due to the complex geomorphic processes driven by Earth's interior and exterior forces. There is still a gap in the numerical modelling of ground landform evolution processes in earth science research.
Here, we report that ruins of an archaeological city in China can act as direct evidence to reconstruct a paleo-surface, which can be further used to facilitate landform evolution modelling. This direct evidence is possible because the discovered archaeological city (named Wucheng, located at 113°33′55″ E, 39°16′31″ N) existed in the Neolithic period and was built on a loess area (Fig. 1). However, during the Han Dynasty (202 B.C.-220 A.D.), this city was abandoned after modern gullies eroded the long sides of the city ramparts 30 . From~2000 years of gully erosion and landform evolution, the current landscape, which has a complete watershed surface morphology, formed. Thus, the archaeological city ramparts are direct evidence indicating the surface morphology of the area, which helps to reconstruct the paleo-surface before the ramparts were eroded. In addition, we simulate the landform evolution process with the reconstructed paleo-surface and the modern surface. Then, the erosion rate is calculated and compared with the climate background information. Finally, the bidirectional interactions between the landform evolution and the rise and demise of human civilisation are further analysed to achieve a better understanding of the human-land relationship.

Results
Reconstructed surface morphology before the Han Dynasty. The morphological characteristics of the archaeological city before the modern erosional gully evolved are quantified with different terrain derivatives ( Fig. 2 and Supplementary Fig. 1). In ancient times, the city contained large smooth areas with low relief surface morphology. Especially in the middle region of this city, no dissected gullies exist that are normally generated by intense erosion. The surface in this area before the Han Dynasty shows much less breakage and fragmentation than the current surface. Because of the landform evolution process that lasted 2000 years, nearly 40% of the city area has been eroded by modern erosional gullies, which now appear as gully-affected areas. The mean elevation of the city decreased by 6.2 m, while that of the gully-affected areas decreased by 12.9 m, with a maximum of 38.4 m. This decreased elevation in the gullyaffected area increased the mean gully-cutting depth from 1.16 m to 2.12 m. The mean slope gradient also increased by 7.22°. All of the above topographical values may even be larger since the city was built upstream of the basin, while the downstream part was likely experiencing more intense gully erosion. In addition, 66,6719 erosional gullies exist in the loess area of the Yellow River basin of China 31 ; if these findings are expanded to the entire loess area, the surface morphology should be much different from that of the present. As the birthplace of Chinese nations, it is difficult Fig. 1 Archaeological ruins of Wucheng city in China. a Location of the archaeological city; it is near the Yellow River and the Great Wall of the Qin-Han Dynasty of China; b detailed surveyed range of the archaeological city and its surrounding geomorphological background; c ancient Wucheng city ramparts; d city rampart constructed by rammed earth; e, f pottery tablets with rope veins from the Han Dynasty discovered in Wucheng. The ramparts were built on a Neolithic settlement lime layer, which can still be found in some places in the northwestern corner of the wall.
to imagine that the current broken and fragmented loess area in the Yellow River basin could have nurtured the Chinese people in ancient times. However, the current study changes the expression and understanding of broken and fragmented surfaces in the current loess area in regard to the pre-Han Dynasty. The pre-Han Dynasty topography was likely more suitable for human survival and development of civilisation in this area. The suitable terrain conditions promoted the development of this ancient city and made the city one of the largest in the region.
Landform evolution process in the last 10 ka. Figure 3 shows the landform evolution process of the archaeological city area during the last 10 ka. According to the simulated surfaces, modern erosional gullies were not developed in the Neolithic age. From ten thousand to five thousand years ago, the whole terrain surface was relatively smooth and had low topographical relief (Fig. 3a, b). With the passage of time, gullies continued to develop and the terrain surface was further divided. As shown in Fig. 3d (where the archaeological city was abandoned during the Han Dynasty), the modelled surface morphology of the watershed was not as broken as that of the modern surface. However, this ancient city had started to be threatened by the development of modern erosional gullies. A modern erosional gully eroded to the foot of the ancient city rampart on the long northwestern side. In addition, incised gullies began to develop on the slopes on both sides of the ancient city (Fig. 3d). The fundamental surface of Wucheng city started to erode. Then, in the next 2000 years, different kinds of gullies developed rapidly, and gully erosion became the main erosion type (Supplementary Figs. 3-5). Figure 3h shows that the statistical indices of the minimum, maximum and average elevations gradually decrease during the landform evolution process. This decreasing phenomenon reveals the dominant erosion force in this area during the last 10 ka. However, the elevation difference of the entire watershed had a significant turning point from a decreasing tendency to an increasing tendency at~2897 B.C. This result suggests that two different dominant erosion forces existed during the landform evolution process of the Wucheng watershed, i.e., hillslope erosion and modern gully erosion. Before 2897 B.C., modern erosional gullies did not reach the watershed, and the entire area did not evolve due to the gully incision process. Then, hillslope erosion was likely the dominant power shaping the surface morphology. With this power, the top hillslope areas with higher elevations were more likely to be eroded, while the bottom areas evolved due to both erosion and accumulation, thus resulting in an overall decreasing tendency for the entire elevation difference (a decrease of 36.3 m).
However, after 2897 B.C., as branches of the Yellow River incised, modern erosional gullies finally reached the Wucheng watershed. The long sides of the city ramparts of Wucheng were eroded, and the city was abandoned. The entire elevation difference of the watershed continued to increase (an increase of 31.4 m to the present), and modern gully erosion was the dominant power shaping the surface morphology in the bottom areas. Because of this dominant power, the modern erosional gully developed rather quickly during this period. Until recently, this quickly increasing elevation difference has stabilised since the modern erosional gully has reached the hard bedrock underlying the loess. The above decreasing and increasing tendencies have also enriched the current understanding of gully erosion, which should be the dominant erosional power in loess areas, transferring loess sediment from upstream to downstream in the Yellow River 32 . This enrichment occurred because, during the above two different tendencies, no remarkable difference in the erosion rate between these two different stages was found (Fig. 4a). In addition, a gully erosion-dominated understanding should rely on whether modern erosional gullies have reached the area or eroded down to the hard bedrock. These findings should also be evaluated by the current landform evolution model of erosional power used in these models, since even in the same concave area within the same watershed, the landform evolution process varies because the erosional power appears different 33 .
Erosion rate curve during the Holocene. Figure 4 shows the variation in the erosion rate of the archaeological city during the Holocene. Overall, an obvious increasing tendency in the erosion rate occurred during the early Holocene (Fig. 4a). The increasing tendency also matched well with the rainfall variability records in this area during this period 34 . In the middle Holocene, this erosion rate tended to fluctuate slightly. When the city was abandoned, the erosion rate tended to increase quickly with a high degree of oscillation, which appears mismatched with the decreasing rainfall variability records. This mismatch was supposedly insufficient human maintenance of the abandoned city, and the natural erosion process would have been accelerated because previous human activities damaged the soil structure in the area. As shown by the detailed erosion rate curve in the last 2 ka (Fig. 4b), most erosion rates range from approximately 3.6 mm/a to 8.2 mm/a, and the average is approximately 6.0 mm/ a. In addition, fluctuations in the speed are exhibited during the landform evolution process. In the last 2 ka, the drivers of the landform evolution process in the loess area were dominated by climate factors due to tectonically stable geological conditions, limited loess dust deposition, a single loess material composition and abandoned human maintenance 35 . Thus, the quick increases and decreases in the erosion rate could be the responses to extreme climatic events in a certain period.
Moreover, as shown in Supplementary Fig. 7, the changing tendency of the erosion rate in different subwatersheds is similar to that of the complete watershed, as shown in Fig. 4b, indicating the inapparent spatial structure of erosion in the Wucheng watershed. However, in the subwatershed located in the abandoned city, the average erosion rate in the last 2 ka appears lower than that in the other watersheds. This result suggests that the remaining anthropogenic ruins still affect the landform evolution process. The area of influence can be large (Supplementary Fig. 7 and Supplementary Table 1), while the erosion depth is limited. In addition, a remarkable decrease in the last 100 years occurs in the curves of all watersheds. During this period, human activities had a remarkable impact on the landform evolution process in this area. Terraces and other anthropogenic facilities have been widely constructed here in the last 100 years 36 , which suggests that current human activities negatively influence and slow erosion on the Loess Plateau.

Discussion
Implications of the landform evolution process for climate change. The results of the landform evolution process in the archaeological city provide the first variation curve of the erosion rate in the loess area during the Holocene, and this erosion rate variation is highly related to climate change. Previous studies have fully investigated paleoclimate change by using materials such as stalagmites or loess strata, and climate change information such as precipitation or temperature during the Holocene has attracted widespread attention among scholars 37,38 . However, the landform evolution process in loess areas should be the consequence of climate change because it is one of the most important drivers. This consequence can also help reveal climate information under the background of stable tectonic conditions, limited loess deposition, and a single material composition in this area. Thus, our erosion rate curve results provide additional evidence indicating climate change during the Holocene. As shown in Fig. 4a, from the last glacial period of the Holocene, the climate of the Loess Plateau tended to be wet and warm, and the erosion rate in this stage began to increase quickly in the early Holocene. Until the Holocene Megathermal, the erosion rate maintained a high degree of oscillation. These oscillations also imply a changing climate background with different landform evolution cycles. This understanding of the climate change effect on the erosion rate can also be supported by recent contributions; e.g., extreme rainstorms in a certain period have been detected as the major driver for the change in erosion rate 39,40 . This means that a remarkable increase or decrease in erosion rates over a certain period usually means that there is a rapid shift between precipitation change environments. In the last 3 ka, especially in the warm period of the Zhou Dynasty 34 , erosion significantly accelerated, and an obvious increase can be found on the erosion rate curve. This finding suggests that the results generated by our model effectively match the climate changes and can support the relationship between loess landform evolution and climate change. This erosion rate curve during the Holocene offers a new perspective for understanding the climate change characteristics of loess areas and even East Asia during the Holocene, which can also help in understanding the environment.
Bidirectional interactions between landform evolution and human civilisation. The reconstructed surface of the archaeological city is the first landscape with clear evidence that indicates the surface morphology before the Han Dynasty. Our results also suggest that this archaeological city was abandoned because of the landform evolution process, which has highly interesting implications regarding interactions between the landform evolution process and the rise and demise of human civilisation. Wucheng is located near the Yellow River and the Great Wall of the Qin-Han Dynasty of China. The conflict between nomadic civilisation and agricultural civilisation continued in ancient China for thousands of years, especially during the Han Dynasty, when the conflict reached a peak 41 . Wucheng is located at the junction of the two civilisations. Wucheng was not only a functioning city but also served the military role of defending against enemies during the Han Dynasty. Thus, the manmade ramparts were quite important for the city during that time. However, with the increased erosion of the gullies, there was a high risk that the city ramparts on the long northwestern side would be destroyed. This condition significantly damaged the military function of the city. In addition, the people of the city faced the high risk of war due to its geographical location, which was another important reason for the abandonment of the city. There are many archaeological sites available around the world, and many of them are caused by geomorphological processes such as erosion and landform evolution [42][43][44][45] . Thus, the idea of using archaeological city ruins as evidence can also be used for landform evolution modelling in other places. For instance, the ancient city Shimao in the loess area was also destroyed by gully erosion 42 , the ancient city Loulan in the desert area of China was badly damaged due to wind erosion and sand accumulation over thousands of years 43 , the ancient city Pompeii in the Campania region of Italy was buried by volcanoes erupting 44 , and historical cities were damaged by avalanches and outburst floods in glacial landform areas 45 . Therefore, at these archaeological sites, the rise and fall of urban civilisation are closely related to the geomorphological process.
In the loess area, smooth topographical conditions with low terrain relief promoted the construction of Wucheng city during the Neolithic period. However, the erosion of the city ramparts in the Han Dynasty caused by intense gully erosion during the landform evolution process made the city no longer suitable for human living; then, the city was abandoned, and the present broken and fragmented surface morphology formed. Previous studies have investigated the impact of human activities on geomorphological processes since the beginning of the Anthropocene 33,46,47 . Such artificial objects should surely affect the Earth's surface processes since they change the original shape of the natural surface morphology. In terms of the impact of the landform evolution process on human civilisation, unlike abrupt geomorphological hazards, the landform evolution process is always slow and occurs over a longer time period. Thus, the impact of the landform evolution process on human civilisation has always been ignored. However, with a better understanding of the past and present surface morphology and by modelling the landform evolution process, we can predict the future surface morphology. Thus, it is necessary to numerically model the landform evolution process of the Earth's surface. Although this process is slow, it can also destroy cities, as shown by this study in ancient times. Ancient Wucheng was once a large and prosperous city in the region. However, with increasing erosion, the developing gullies severely damaged the surface of ancient Wucheng and finally led to the failure and abandonment of ancient Wucheng. Many ancient cities exist in loess areas, and many of them might have faced a historical situation similar to that of ancient Wucheng and specifically were threatened by intense erosion and the development of gullies. People who lived in these cities had to abandon their homes and migrate.

Methods
Data acquisition. The location of the archaeological city of Wucheng (39°16′49.2″ N, 111°33′ 49.3″E) is in Shanxi Province, China, which exhibits large areas of typical loess landforms. This nearly rectangular archaeological city is also near both the Yellow River and the Great Wall of the Qin-Han Dynasty of China. The archaeological city of Wucheng is one of the Key Cultural Relics Protection Units in Shanxi Province 30 . The rampart length in ancient Wucheng is~1000 m in the long direction and~500 m in the short direction. The northern and southern ramparts were built on the top of a hill along ridges. The eastern and western ramparts were built gradually down the slope to the bottom. According to the discovered relics and antiques, the survival period of ancient Wucheng was approximately from the Neolithic period (~10 ka B.P.-5 ka B.P.) to the Han Dynasty (202 B.C.-220 A.D.). Due to the landform evolution process, modern erosional gullies damaged ancient Wucheng and caused the destruction of long sides of the city ramparts, which led to the abandonment of ancient Wucheng in the Han Dynasty. The ruins of the short side of the city ramparts still exist and act as the watershed boundary.
Since the archaeological city was abandoned in the Han Dynasty, the modern erosional gully had not broken the city ramparts before the Han Dynasty. Thus, the paleo-surface before the modern erosional gully evolved can be reconstructed with the preserved city information that acts as the archaeological ruins of the city. The archaeological ruins of the city used in this study contain samplings along preserved ramparts and hillslope areas that were not eroded by modern erosional gullies. The basic elevation data for landform evolution modelling were collected by unmanned aerial vehicle (UAV) photogrammetry technology. Stereo pairs shot by the camera carried by a UAV were transformed into point cloud data and then used to generate a raster digital elevation model (DEM) 48,49 . In addition, we used real-time kinematic global positioning system (RTK-GPS) devices to collect the location and elevation of sampling points located near the foot of the ancient ramparts and relics. Then, spline interpolation was applied to reconstruct the fundamental surface of ancient Wucheng on the basis of sampling points. Moreover, archaeological data, including historical literature and pottery clay pieces, were collected through field investigations. These data can help define the existence and abandonment period of ancient Wucheng and further construct a time-series DEM sequence.
Landform evolution model. The landform evolution model is a physical simulation of geomorphological processes via computational techniques 49,50 . Many factors influence the landform evolution process, and the general form of the conservation equation of the landform evolution process can be expressed as follows: where U is the rock uplift rate; K is the erosion coefficient; A is the upstream drainage area; S is the slope factor; m and n are the exponents; C is the erosion threshold; D is the coefficient of diffusion; z is the elevation; and t is the time.
In the last 10 ka, the geological background of the archaeological city appears tectonically stable with a single material composition, namely, loess. Thus, rock uplift has limited influence on landform evolution on the Loess Plateau. In addition, the vegetation of this area has been relatively sparse in the past 2000 years due to the climate conditions 51,52 , which further intensifies soil erosion. Thus, the erosion process dominated the landform evolution process during this period. Therefore, it can be concluded that the advancement of gullies under the action of erosion is an important reason for the destruction of Wucheng. In other words, factors A and S are the key factors that should be considered in the evolution model, which is also in accordance with the soil erosion model of the Revised Universal Soil Loss Equation (RUSLE) in the loess area. In addition, the coefficients and exponents in the model are difficult to measure directly 53,54 . To avoid this difficulty, we integrated the intelligence method to simulate the landform evolution process with the key factors of A and S.
The intelligence method is a spatiotemporal dynamic simulation system (SDSS) that consists of four elements: cells, states, neighbours and rules 55,56 . Cells are basic units of cellular automata, while states are the attribute values of cells at a certain time. These cells that interact and influence each other gradually transform into new states and simulate the evolution of landforms as the SDSS method runs. In addition, in the SDSS-based landform evolution model, we established the cell transition rule on the basis of the key factors of A and S, and these two key factors were modified according to the RUSLE. This is because the RUSLE has been suggested to be more suitable for loess areas when estimating sheet erosion and gully erosion 48,57 . The RUSLE is a classic erosion model, and its effectiveness in the Loess Plateau has been extensively explored and validated 58 . The calculation of the S factor concerns the terrain slope gradient and soil type and is generally different in different regions. According to the erosion type and soil properties in the Loess Plateau, several empirical formulae have been proposed to express the relationship between the terrain slope gradient and the S factor. We chose the most popular Eq.
Then, factor A was modified as factor L to adapt to the erosion process in the Loess Plateau. Factor L defines the impact of slope length and controls the degree of soil erosion 58 . The calculation of factor L can be followed by Eq. (4): where A is the upstream drainage area, R is the cell size, and m is the coefficient of slope gradient and can be obtained by Eq. where β is the terrain slope gradient. Thus, by considering the key factors of A and S in the landform evolution process, the SDSS-based landform evolution model adjusts the elevation of each pixel to generate the terrain surface of the next period. The workflow and cell status in the SDSS-based landform evolution model can be found in Supplementary  Figure 8 and Supplementary Table 2, respectively. Finally, we designed a time calibration rule to construct the connection between the number of iterations and the historical time and determine the relatively correct period of DEMs. The correlation coefficient ρ was calculated to quantify the similarity between the fundamental surface of ancient Wucheng S 1 and the reconstructed result S 2 . Cov where CovðS 1 ; S 2 Þ is the covariance of two surfaces; D is the variance; and E is the mathematical expectation. The value domain of ρ is [−1,1], and the higher ρ is, the higher the similarity between two surfaces. Then, we calculated several indices that can indicate the surface characteristics and analyse the landform evolution in the past 10 ka.
where D is the cutting depth, E ave is the average elevation in a specific analysis window, and E min is the minimum elevation in the same area.

Rough ¼ 1 cosðβÞ
where Rough is roughness, and β is the slope gradient of each pixel.
where ER is the erosion rate; E t ave and E tÀ1 ave are the average elevations of the DEM of periods t and t-1, respectively; and 4t is the interval between these two periods.

Data availability
The data related to this paper can be found on the Figshare website (https://figshare.com/ articles/journal_contribution/landform_evolution_modelling/21581220).

Code availability
The codes related to this paper can be found on the Figshare website (https://figshare. com/articles/journal_contribution/landform_evolution_modelling/21581220). The information of model parameters can be found in explanation file in this repository.