WANG Yao, HOU Li-sheng, CAI Yun-long

1 Development Research Center of China Geological Survey, Beijing 100037, China.

2 China Geological Survey, Beijing 100037, China.

3 Center for Land Study, College of Urban and Environmental Sciences, Peking University, Beijing 100871, China.

Abstract: In recent years, research on spatial scale and scale transformation of eroded sediment transport has become a forefront field in current soil erosion research, but there are very few studies on the scale effect problem in Karst regions of China. Here we quantitatively extracted five main factors influencing soil erosion, namely rainfall erosivity, soil erodibility, vegetative cover and management, soil and water conservation, and slope length and steepness. Regression relations were built between these factors and also the sediment transport modulus and drainage area, so as to initially analyze and discuss scale effects on sediment transport in the Wujiang River Basin (WRB). The size and extent of soil erosion influencing factors in the WRB were gauged from: Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM), precipitation data, land use, soil type and Normalized Difference Vegetation Index (NDVI) data from Global Inventory Modeling and Mapping Studies (GIMMS) or Advanced Very High Resolution Radiometer (AVHRR), and observed data from hydrometric stations. We find that scaling effects exist between the sediment transport modulus and the drainage area. Scaling effects are expressed after logarithmic transformation by a quadratic function regression relationship where the sediment transport modulus increases before decreasing, alongside changes in the drainage area. Among the five factors influencing soil erosion, slope length and steepness increases first and then decreases, alongside changes in the drainage area, and are the main factors determining the relationship between sediment transport modulus and drainage area. To eliminate the influence of scale effects on our results, we mapped the sediment yield modulus of the entire WRB, adopting a 1 000 km2 standard area with a smaller fitting error for all sub-basins, and using the common Kriging interpolation method.

Keywords: Sediment transport modulus; Scale effect; Soil erosion; Wujiang River Basin

The phenomenon of eroded sediment transport occurs in a certain time and space range, and is the result of mutual action and restriction of various natural and social factors. In drainage basins with different scales and different areas, factors (landform, precipitation, vegetation, soil, land use, etc.) that influence eroded sediment transport are complex and varied. In recent years, many scholars have noticed the change in spatial scale in soil erosion research, and research on spatial scale and scale transformation of eroded sediment transport has become a forefront field (NI Jiu-pai et al. 2010; YAN Yun-xia et al. 2010; YUAN Zai-jian et al. 2007; LIAO Yi-shan et al. 2008; Kirkby M J and Mcmahon M L, 1999; Kirkby M J and Cox N J, 1995; Kirkby M J et al. 1996; Favis-Mortlock D T et al. 1996; Renard K G et al. 1997) in current soil erosion research. Based on the study on relations between drainage area and three variants, namely sediment amount, erosion amount and sediment transport ratio of the Ganjiang River Basin, JING Ke et al. (2010) found no significant correlation between the three variants and drainage area, and factors influencing the three variants had no scale effects on drainage area. Research on the upper reaches of Yangtze River found that changes in rainfall erosivity and soil erodibility factors along with changes in drainage area, caused a negative power function unit and regression relationship between the sediment transport modulus and the drainage area (SHI Chang-xing, 2008). However, clear differences in the relationship between sediment transport modulus change and the basin scale of the main branches the upper Yangtze River, indicate many unresolved issues requiring further in-depth study. Karst regions form the basis for much of the focused research on soil erosion in China to date. Furthermore, much research has centred around soil erosion, sediment transport and yield (XIONG Kang-ning et al. 2012; WAN Jun et al. 2004; XU Yue-qing et al. 2006; WANG Wen-bo, 2008; ZENG Ling-yun, 2008; XU Yue-qing et al. 2008; WANG Hong-ya et al. 2006; ZHU An-guo et al. 1993; CHEN Song-sheng et al. 2008; XIONG Ya-lan et al. 2011; HE Yong-bin et al. 2009) but there are very few studies on the scale effect problem. This paper analyzes and discusses the scale effect on eroded sediment transport in Wujiang, a main branch of the upper Yangtze River located in southwest Yunnan- Guizhou Plateau, including the influence of karst landforms. This work is highly significant in illustrating eroded sediment transport mechanisms and establishing a regional soil erosion model.

Fig. 1 River and elevation map of Wujiang River Basin in Guizhou Province, in southwest China

1 Research area overview

Our research area is the Wujiang River Basin (WRB) (Fig. 1), running through the middle of Guizhou Province. Wujiang River is a branch of the Changjiang River system, originating from Xianglu Mountain in Weining County at the eastern foot of Wumeng Mountains to the west of Guizhou Province. For the majority of its' length, the WRB flows through 23 counties and cities such as Xianning, Shuicheng and Nayong, and enters Sichuan Province from the northeast of Guizhou Province. The main branches of the Wujiang River include the Liuchong River, the Sancha River, the Maotiao River, the Yachi River, the Xiang River, the Pianyan River, the Qingshui River, the Shiqian River and the Yuqing River. Flowing from the west to east, the WRB covers 51 817 km2in Guizhou Province, of which 75.6% is area characterised by karst landforms developed by carbonate rocks. As they are not affected by headward erosion, upstream areas Weining and Hezhang maintain good plateau ground, while middle and downstream areas are characterised by a wealth of deep valleys and high mountains, with surface relief reaching as high as 300 m.

2 Research data

A wide range of data was collated for the WRB study area, including the Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER GDEM) data, 1/1 000 000 soil type data, GIMMS/AVHRR NDVI data during the period 1990-2000, and daily precipitation data from 18 meteorological stations in the WRB and surrounding areas between 1990 and 2000. ASTER GDEM data were downloaded from the Earth Remote Sensing Data Analysis Center, soil type data were downloaded from Environmental and Ecological Science Data Center for West China, daily precipitation data from meteorological stations were downloaded from China Meteorological Scientific Data Shared Service Network, GIMMS/AVHRR NDVI data were from the GIMMS/NDVI data set published by NASA's Global Inventory Modeling and Mapping Studies (GIMMS) and downloaded from Environmental and Ecological Science Data Center for West China (temporal resolution was 15 d, and spatial resolution was 8 km×8 km).

3 Research method

Factors influencing soil erosion were primarily calculated using ArcInfo, Excel and Microsoft Access with reference to the Revised Universal Soil Loss Equation (RUSLE) (Wischmeier W H and Smith D D, 1987).

3.1 Rainfall erosivity factor R

The experimental result of soil loss on karst slopes indicated that surface runoff was generated when rainfall reached 10.4 mm and when the site had experienced earlier sustained rainfall. Clear runoff appeared when precipitation was greater than 15 mm and there had been no earlier rain. Therefore, this paper adopts 12 mm as the critical erosive rainfall value. The calculation of rainfall erosivity required adopting a simple model of monthly precipitation erosivity built by YU B et al. (1998) (Vol N, 1996; Yu B, 1998) in Australia on the basis of daily precipitation data:

where Ejis the precipitation erosivity of month j, MJ·mm·hm-2·h-1; Rdis the daily precipitation mm, and R0is the erosive daily precipitation intensity threshold value, f -frequency, taking 1/12; ω - 5π/6; N - the number of days with rainfall surpassing R0in the j month; S- rainfall during May-August in the summer half year; P- annual average rainfall mm; and α, β and η were model parameters. The parameter β ranges from 1.2 to 1.8, but 1.5 was adopted in this paper according to previous research result and combined the practical condition of the WRB. The parameters α and β met in Equations 1 and 2 when annual rainfall was greater than 1 050 mm, and they met in equations 1 and 3 where annual rainfall was greater than 500 mm but smaller than 1 050 mm. The parameter η was calculated according to annual average rainfall (Equation 4).

3.2 Soil erodibility factor K

Soil erodibility refers to soil sensitivity to erosion, or complexity of soil being separated, washed and transported by rainfall erosion. Acquiring an actual measured K value for different soil types in the WRB is helpful to accurately estimate the amount of soil erosion. However, as there is no condition to acquire an actual measured K value, the K value has to be determined with the assistance of relationships built around the soil erodibility and soil property parameters. This paper estimates K value using factors involving soil organic matters and particles with the help of a widely used set of equations adopted by Williams et al. (Williams J R et al. 1983) in the EPIC model, as follows:

In these equations, K is soil erosion resistance, for which the British unit is t·ac·h/(100 ac·ft·t·ih), and the metric unit is m·t·hm2·h/(hm2·MJ·mm). The metric K value is equivalent to 0.13* of the British K value. Data adopted in this paper are from the Harmonized World Soil Database (HWSD) established by FAO and Vienna International Institute for Applied Systems Analysis (IIASA). For our purposes, the soil classification system largely adopted that used by FAO-90, and American classification standards were adopted for soil texture classification. msrepresents sand content (0.05-2.00 mm), msiltrepresents silt content (0.002-0.005 mm), mcrepresents clay content (<0.002 mm), and orgCρ represents organic carbon content.

3.3 Vegetation cover and the manage ment factor C

Vegetation cover and the management factor in RUSLE is defined as the ratio of loss in specific vegetation-covered land and the loss of tilled clean fallow land under the same conditions. The normalized NDVI index appropriately reflects vegetation coverage. In the image, the NDVI value of each pixel element can be viewed as a weighted average of the NDVI of an area with vegetation cover and the area without vegetation cover. Vegetation coverage per pixel element is represented by the weight of the NDVI value of the area with vegetation cover, while the NDVI value of the area without vegetation cover represents the “I-vegetation coverage” (CHENG Sheng-dong, 2010). The extent of vegetation coverage in the image can be calculated using the following equation:

In the equation, Vnrepresents vegetation coverage of the pixel element; NDVImaxstands for the NDVI value of bare soil or area without vegetation coverage (in other words, the NDVI value of the pixel element without vegetation); and the NDVIminrepresents the NDVI value of the pixel element covered entirely by vegetation.

The value of C is then calculated using the regression equation of the relationship between the C factor value and the vegetation coverage established by CAI Chong-fa et al (2000).

The Maximum-Value Composite (MVC) was obtained for the half-month NDVI data from 1990-2000 adopted in this paper (Holben B N, 1986). The MVC can further eliminate partial interference from cloud, atmosphere and solar altitude. The principle is selecting maximum reflectivity of each pixel element at the composite phase, to reflect maximum photosynthesis of vegetation of each pixel element at the composite phase. Following data smoothing, the half-month NDVI data were subjected to the monthly maximum composite and the annual maximum composite respectively, from which the monthly maximal NDVI (MMNDVI) and the annual maximal NDVI (AMNDVI) were obtained. The NDVI of each pixel element was subjected to maximum composite calculation according to the following equation.

The MNDVIiis the maximal NDVI of the i month obtained by compositing the two half-month images of the i month, and is used for reflecting vegetation cover conditions of each pixel element in this month under the best weather condition. The AMNDVI is calculated using the maximum compositing of MNDVIiover 12 months of a certain year, reflecting the NDVI value of each pixel with a high density of vegetation.

3.4 Soil and water conservation factor P

Without plot experimental data for soil and water conservation factors in the WRB, we estimated soil and water conservation factors using empirical values. Instead, this paper assigns land-use type in the WRB with corresponding soil and water conservation factor values defined by field investigations in the WRB, and combined related research results of XU Yue-qing et al. (2006) in the Maotiao River Basin.

Table 1 P values of different land-use types, namely, land cover type, water unused land, land under urban construction, paddy field, dry land, woodland, grassland

3.5 Slope length and steepness factor LS

Increasing steepness is paired with a reduction in soil particle slide resistance and an increase in the kinetic energy of flows. Accordingly soil loss amount increases. Furthermore, the longer the slope length, the more substantial the total amount of soil loss. Slope length and the steepness factor LS refer to the ratio of soil loss on slopes with a certain steepness and slope length. Soil loss of a typical slope with a standard runoff plot with some other conditions. The factor LS is the accelerated factor of erosion force. With rapid development of GID technology, the pixel element based slope length and steepness calculation have gradually replaced the pattern spot based slope length and steepness calculation. Slope length refers to maximum ground distance from one point on the ground (one grid unit in the DEM) to the starting point of the flow paths. When we combine the data of the flow starting point and flow paths, the maximum cumulative flow length between each grid unit to the starting point along the flow paths is the slope length from the grid to the slope crest.

Both slope length L and steepness S are calculated according to the equation defined McCool et al (Renard K G et al. 1997):

In the equation, L is the slope length factor; 1 is the slope length of the pixel element in feet (ft); 72.6 is the slope length of the American standard runoff plot (ft); and m refers to the slope length index, the values of which are as follows:

Table 2 Values of slope length exponent (m) (McCool et al. 1997 (Renard K G et al. 1997))

The calculation equation for the steepness factor is as follows:

Calculation of the slope length factor on soil erosion mechanisms is complicated, and there is little in the literature to assist. The ArcInfo-based AML calculation script provided by Robert J. Hickey from University of Washington was downloaded from the website and run using Windows-based ArcInfo Workstation. The work's root catalog, complete work path, grid dem name, basin boundary grid layer name, DEM unit and steepness segmental values were entered as per the system prompts. As steepness was entered in degrees, the value range of the segmental value fell between 0-1 degrees. The smaller the segmental value, the more accurate the slope length calculation, but the level of calculation increased. Hence, it is recommended that the segmental value is set at 0.7 degrees when steepness is smaller than 5 degrees, and should be set at 0.5 degree when steepness is larger than 5 degrees. The final generated value slp_lgth_ft is the cumulative flow length, namely slope length (in feet). The value c_ruslels1 is the slope length L factor, c_ruslels is the steepness factor, and c_ruslels2 is the slope length and steepness factor LS.

4 Results and analysis

4.1 Analysis of the scale effect on sediment transport in the basin

The plotted relationship between sediment transport modulus and drainage area recorded in the control zones of the 11 hydrological stations in the WRB is shown in Fig. 2. In the log-log coordinate plot, the x coordinate is the drainage area S controlled within the stations and the y coordinate is the sediment transport modulus A. From the regression relationship analysis, we find that there is a significant correlation in the quadratic function relationship of the drainage area and the sediment transport modulus when the quadratic function relationship is 0.63. In formula, this relationship is represented in the following equation, where S refers to the sediment transport modulus (t·km-2·a-1), A is the drainage area (km2), r2is the correlation coefficient, and n represents the sample number:

Fig. 2 Relationship between sediment transport moduli and drainage areas recorded at hydrological stations in the Wujiang River Basin

Fig. 2 illustrates a sediment transport scale effect, whereby the sediment transport modulus of the basin shows an increasing trend before decreasing alongside increases in the drainage area. To explore the in-depth reasons for this, the comprehensive relationships between weather, landform, and surface material composition and that between R, C, K, LS, and P of land use factors are considered in relation to changes in the drainage area. It can be seen from Fig. 2 that the basin's sediment transport modulus and drainage area are correlated, tending to increase before decreasing, and scale effects exist on sediment transport.

4.2 Analysis of scale effects on the influencing factors of soil erosion

Arc Hydro Tools generated 189 sub-basin samples (Fig. 3) within the 150 km2catchment area. The statistical relationship between the catchment area and abovementioned five soil erosion influencing factors was calculated, and scale effects of each influencing factor was analyzed.

A spatial distribution map is presented below, illustrating the annual average rainfall erosivity factor R, soil erodibility factor K, slope length and steepness factor LS, vegetation cover and management factor C, and soil and water conservation factor P, representing the climate, landform, surface material composition and land use factors from 1990 to 2000. A table and map demonstrating the relationship between R, K, C, P and LS along with change of drainage area are also presented here.

Fig. 3 Spatial distribution maps of the Revised Universal Soil Loss Equation Factors P, R, K, LS and C in the Wujiang River sub-basins, China

Table 4 Regression relationships between each Revised Universal Soil Loss Equation factors and drainage area in the Wujiang River Basin, China

Fig. 4 Regression relationships between the Revised Universal Soil Loss Equation factors and the drainage area

Among the five factors influencing sediment transport, the terrain factor LS shows an increasing trend before decreasing, associated with an increase in drainage area. Soil erodibility factor K shows an increasing trend, vegetation cover and management factor C shows a decreasing trend before increasing, soil and water conservation factor shows no obvious change; and annual average rainfall erosivity factor R shows a decreasing trend. Therefore, the deduction is that the change in the terrain factor LS leads to change in sediment yield modulus in the WRB, which increases first before decreasing as the drainage area increases. To do this, the sediment yield moduli of all basin intervals were converted to moduli under a certain standard area, and a sediment yield modulus interpolation calculation was performed. Referring to research methods of others in the field, this paper attempted a specific scale transformation method by combining the practical condition of the WRB, and thereby formulated a scale effect based sediment yield modulus distribution map.

4.3 Correction of scale effects

Drawing a sediment transport modulus map traditionally involves selecting several hydrometric stations where a long-series of sediment yield modulus observation has occurred, determining the drainage basins or the basin intervals represented by the hydrometric stations, and then drawing the sediment yield modulus isoline using an interpolation method. Given that such a sediment yield modulus map involves different basin areas which may have considerable influences on the sediment yield modulus, scale effects on the basin areas from the sediment yield modulus must analyzed in order to perform a scale transformation.

4.3.1 Correction method

From the above analysis, it can be seen that scale effect calculating the WRB's sediment transport modulus involves consideration of scale effect. What this correction means is that is that the sediment yield modulus will be displayed at the same scale on the distribution map (for example 1 000 km2) to eliminate the influence of the scale effect. Suppose the fitted equation of the quadratic equation Y=AX2+BX+C. If the practical data of a certain station to be corrected is (x1, y1), and the correction equation is Y'=A(X')2+BX'+C', then C'=y1-(AX12+BX1), which can be converted to Y'=A(X')2+B X'+(y1-AX12+BX1), and thenWhen the drainage area Adof the drainage basins controlled by the hydrometric stations have different values (such as 1, 100, or 10 000 km2), the corresponding corrected sediment yield modulus can be obtained.

4.3.2 Drawing the sediment transport modulus map

According to the correction method, each sub-region's sediment transport modulus was corrected as per the standard area, and the research area's sediment transport modulus map was drawn using the spatial interpolation method. When the data demonstrated correlations within a certain distance or in a certain direction, and could reflect a certain trend in the changing space, the Kriging spatial interpolation method was the best choice. The range of the station area included in the sub-regions and the fitted equation were different, resulting in different optimum areas. Using the common Kriging interpolation method, our work drafted a sediment transport modulus map for the whole research area by adopting a 1 000 km2scale, with a smaller fitting error for all sub-basins as the standard area.

Fig. 5 The average annual sediment transport modulus map corrected to the standard area of 1 000 km2 in the Wujiang River Basin, China

The distribution of the annual average sediment transport modulus at the erosion interval of 230-450 t/(km2·a) did not change greatly in the study area, and was especially consistent before and after the scale correction (Fig. 5 and Fig. 6). However, the upstream and downstream of the western WRB changed significantly pre and post correction, as did the average annual sediment transport modulus of the Sancha River in the west and the upstream area of the Liuchong River, and in the WRB's downstream area behind the Jiangjie River.

Fig. 6 The average annual sediment transport modulus map of the Wujiang River Basin, China

5 Conclusions

(1) Scale effects exist in the relationship between the WRB sediment yield modulus and the drainage area. The sediment transport modulus and the drainage area form a quadratic regression relationship after logarithmic transformation, demonstrating that the sediment transport modulus increases before decreasing, alongside changes in drainage area. Among eroded sediment transport factors, the effect of slope length and steepness both increase and then decrease, in relation to changes in drainage area. This result indicates that slope length and steepness are the major factors determining the relationship between the sediment transport modulus and the drainage area.

(2) The drainage sediment yield modulus map is drawn as per the standard area of 1 000 km2after scale correction. This map illustrates is remarkable difference in the annual average sediment transport modulus in the upstream areas of the Liuchong and Sancha Rivers, as well as Wujiang River's downstream area behind Jiangjie River. In particular, after correction the sediment transport modulus of the upstream areas decreases, while that of downstream area increases. Before and after correction, the sediment transport modulus' high value center is concentrated in the upper Wujiang River in the Yachi River drainage system, including Dafang County, Nayong County, Zhenxiong County and Zhijin County. What this means is that this sediment yield region urgently requires treatment.

(3) In the past, purely considering the relationship between the sediment transport modulus and drainage area or the influence of individual factors does not give us a full picture of soil erosion mechanism. We have instead made a comprehensive comparison and quantitative analysis of the relationship between several key impacting factors of soil erosion and the drainage area, which reveal in depth how scale effects on sediment transport in the WRB are generated. Our research result aids improved understanding of the general law of sediment transport in drainage areas, as well as treatment options for water and soil loss. However, according to our scaling effect analysis on more than 180 sub-basins in the WRB, the changing scope of the drainage area ranges from several square meters to around 1 000 square meters, a result which differs from the change in scope of the drainage area obtained through measured data during the scaling effect analysis. There could be a variety of reasons for this difference, all requiring further in-depth study. In addition, along with consecutive large-scale water and soil loss treatment works carried out in the WRB, water conservancy projects have developed continuously, and vegetation cover, land use, soil and water conservation and other factors may change remarkably. The sediment yield and transport law of the WRB may therefore change in correspondence, thereby requiring further monitoring and research with the support of continuously updated materials.

Acknowledgements

This study is generously supported by Project of National Natural Science Foundation of China (41641011), National Geology and Mineral Resources Survey and Assessment Program (DDT0160087).