Emergent stationarity in Yellow River sediment transport and the underlying shift of dominance: from streamflow to vegetation

Soil erosion and sediment transport play important roles in terrestrial landscape evolution and biogeochemical cycles of nutrients and contaminants. Although discharge is considered to be a controlling factor in sediment transport, its correlation with sediment concentration varies across the Yellow River basin (YRB) and is not fully understood. This paper provides analysis from gauges across the YRB covering a range of climates, topographic characteristics, and degrees of human intervention. Our results show that discharge control on sediment transport is dampened at gauges with large mean annual discharge, where sediment concentration becomes more and more stable. This emergent stationarity can be attributed to vegetation resistance. Our analysis shows that sediment concentration follows a bell shape with vegetation index (normalized difference vegetation index, NDVI) at an annual scale despite heterogeneity in climate and landscape. We obtain the counterintuitive result that, as mean annual discharge increases, the dominant control on sediment transport shifts from streamflow erosion to vegetation retardation in the YRB.

1 Introduction

Watershed sediment transport, from hillslope to channel and subsequently the coast, is crucial to erosion management; flood control; river delta development; and the quantification of global biogeochemical cycles of materials such as organic phosphorus, iron, and aluminum (Martin and Meybeck, 1979; He et al., 2014). During the 20th century, human activities significantly modified the landscape, leading to a reduction in sediment yield and coastal retreat worldwide (Walling and Fang, 2003; Syvitski et al., 2005). Known for its severe sediment problems, the Yellow River (YR) has been a hotspot for studies on soil erosion and sediment transport for decades. Since the 1950s, the annual sediment yield has declined by 80 % because of check-dam construction and ecosystem restoration such as the Grain-for-Green project, motivating discussion on the necessity for further expansion of revegetation schemes (Chen et al., 2015).

Most studies on the physical mechanisms of soil erosion and sediment transport were conducted in relatively small sub-catchments (Collins et al., 2004; Ran et al., 2012). In order to interpret the patterns discovered at the basin scale, it is essential to understand the scaling effects of soil erosion and sediment transport. Specifically, would the mechanisms identified at a small scale also prevail at the basin scale? If not, what factors influence upscaling (Mutema et al., 2015; Song et al., 2016)? However, existing studies on the scaling effects of sediment transport are rather limited and show no significant spatial coherence in the scaling of sediment transport (Le Bissonnais et al., 1998; Deasy et al., 2011; Song et al., 2016). Due to the great heterogeneity in the YRB, scaling patterns could be different even within one tributary. Taking the Wuding River as an example, event mean concentration could decrease downstream after the initial increase in one sub-catchment (Zheng et al., 2011) or keep rising until reaching a plateau in another sub-catchment nearby (Fang et al., 2008). Not only the sediment concentration, but also its correlation with discharge varies across the YRB. Although discharge is considered to be one of the controlling factors in sediment transport, how its influence upscales remains to be fully understood. Therefore it is necessary to expand our findings concerning sediment transport from single tributaries to larger scales, especially incorporating diverse climate, environmental, and anthropogenic characteristics, so that we can derive an understanding applicable to the whole YRB. In this paper, we collected observations across the Yellow River basin (YRB) to quantify changes in sediment concentration in the recent decades (Rustomji et al., 2008; Miao et al., 2011; Wang et al., 2016). By analyzing data from gauges across the YRB (Fig. S1 in the Supplement), we attempt to understand how the correlation between sediment concentration and discharge varies across spatial and temporal scales, what are the dominant factors influencing sediment transport in the YRB, and how their contributions vary from place to place.

2 Data and methodology

We collected daily discharge and sediment concentration data from 123 hydrology gauges within our study area: the YRB above the Sanmenxia station, which is the major hydropower station on the YR. From these we selected 68 gauges spanning a range of climate conditions and physiographic areas, from the gauge at the most upstream end of the main stem to the gauges above Tongguan, which is just 100 km upstream of Sanmenxia Dam (Fig. S1). These gauges were selected for at least 15-year (1971–1986) continuous daily discharge and sediment concentration records between 1951 and 1986. For comparison and further examination of our hypothesis, we also extract the annual discharge and concentration data between 2000 and 2012 for seven gauges located at the outlet of the major tributaries from the Yellow River Sediment Bulletin (Fig. S1 green stars).

The vegetation data used in this study correspond to the normalized difference vegetation index (NDVI), which is an index calculated from remote sensing measurements to indicate the density of plant growth (Running et al., 2004). The NDVI data were downloaded from NASA's Land Long Term Data Record (LTDR) project, which provides daily NDVI observations globally at a spatial resolution of 0.05^∘. Instead of the NDVI obtained from Global Inventory Modeling and Mapping Studies (GIMMS), the NDVI from the LTDR project is chosen for its better estimation in the YRB (Sun et al., 2015). The daily NDVI data from 44 gauges located on the eight major tributaries were collected and extracted according to the drainage area of the study gauges from 1982 to 2012 (Fig. S1 green stars). The gauges on the main stem of the YR were not used as the water and sediment condition there is more likely controlled by the major dams along the main stem rather than the hillslope characteristics. Annual maximum NDVI values were used to represent the highest vegetation productivity. The precipitation and leaf area index (LAI) data of the US catchments used for comparison are assembled from the first author's previous work (Ye et al., 2015).

To examine the coupling between discharge and sediment concentration at various temporal scales, wavelet coherence analysis was applied to the daily discharge (m³ s⁻¹) and sediment concentration (kg m⁻³) data following Grinsted et al. (2004). Wavelet transforms decompose time series into time and frequency and can be used to analyze different parts of the time series by varying the window size. They have been applied to geophysical records for the understanding of variability at temporal scales. To examine the co-variation between discharge and concentration in the time frequency domain, we used a wavelet coherence defined as (Grinsted et al., 2004)

$$\begin{array}{}\text{(1)}& R^{\mathrm{2}}\left(s\right)={\displaystyle \frac{\mathrm{|}S\left(s^{-\mathrm{1}}{W}^{XY}\left(s\right)\right){\mathrm{|}}^{\mathrm{2}}}{S\left(s^{-\mathrm{1}}{\left|W^X\left(s\right)\right|}^{\mathrm{2}}\right)\cdot S\left(s^{-\mathrm{1}}\mathrm{|}{W}^Y\left(s\right){\mathrm{|}}^{\mathrm{2}}\right)}},\end{array}$$

where S is a smoothing operator, W^XY is the cross-wavelet transform of time series X and Y representing the common power between the two series, s refers to the scale, and W^X and W^Y are the continuous wavelet transforms of time series X and Y respectively. The wavelet coherence can be considered a correlation coefficient of the two time series in the time frequency domain. The region of cone of influence (COI) was delineated in the wavelet coherence images to avoid reduction in confidence caused by edge effects. Localized wavelets were also averaged through temporal scales to obtain global wavelet coherence (Guan et al., 2011). More detailed explanation about wavelet coherence analysis can be found in Grinsted et al. (2004).

The annual discharge (Q_a) and the sediment yield (L_a) were aggregated daily to further examine their correlation:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{a}}=\left({\sum }_{i=\mathrm{1}}^n(Q_i\cdot \mathrm{3600}\cdot \mathrm{24})\right)/A_{\mathrm{d}}\cdot \mathrm{1000},\text{(3)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{a}}=\left({\sum }_{i=\mathrm{1}}^n(Q_i\cdot C_i\cdot \mathrm{3600}\cdot \mathrm{24})\right),\end{array}$$

where Q_i (m³ s⁻¹) and C_i (kg m³) are the daily discharge and sediment concentration, A_d is the drainage area (km²) of each gauge, and n is the number of days in each year. This analysis is applied only at an annual scale since this is when the coupling from wavelet coherence analysis is the strongest (the one with the largest wavelet coherence). The annual mean concentration (C_a) was calculated as

$$\begin{array}{}\text{(4)}& C_{\mathrm{a}}=L_{\mathrm{a}}/(Q_{\mathrm{a}}\cdot A_{\mathrm{d}}/\mathrm{1000}).\end{array}$$

The long-term mean annual discharge (Q_m) and the long-term mean annual concentration (C_m) were also calculated by averaging for the period of 1951 to 1986. Note that both the parameters Q_a and Q_m used here are area-specific discharges (mm yr⁻¹). For each gauge, a linear regression was fit to describe the correlation between annual discharge (Q_a) and annual mean concentration (C_a). The slope of this linear regression (α_QC) is used to describe the rate of change in sediment concentration with changing discharge at an annual scale.

Figure 1 Scatter plots between long-term mean annual discharge (Q_m) and (a) wavelet Q–C coherence at daily, monthly, and annual scales from the 68 study gauges and (b) slope of the discharge–sediment-concentration regression (α_QC) at an annual scale from the 68 study gauges (R²=0.55 and p value <0.0001).

3 The emergent stationarity in sediment concentration

We applied wavelet coherence analysis to daily discharge and sediment concentration data at 68 study gauges across the YRB (Figs. S2, S3). The results show that, across the gauges, the coupling between discharge and concentration (Q–C) declines with mean annual discharge (Q_m) at all three temporal scales (Fig. 1a). That is, as Q_m increases, the influence of streamflow on sediment transport becomes weaker and weaker across the gauges, both at intra-annual and within-year scales.

This fading impact of streamflow as it increases can be further quantified in terms of a linear regression between discharge (Q_a) and mean sediment concentration (C_a) at an annual scale, when the coupling between discharge and concentration (Q–C) is the strongest (Fig. S4). As can be seen from Fig. 1b, though annual mean concentration is positively correlated with annual discharge at most gauges, the slope in the Q–C regression (α_QC) declines exponentially with Q_m across the gauges (p value <0.0001). The larger Q_m is, the less sensitive sediment concentration responds to variation in annual discharge. For example, gauges with α_QC less than 0.1 are the ones with Q_m larger than 60 mm yr⁻¹. When Q_m is larger than 100 mm yr⁻¹, the variation in sediment concentration is less than 1 % of that in streamflow (α_QC<0.01), and thus sediment concentration can be approximated as invariant to changing discharge. Most of these gauges are located on the main stem or near the outlets of tributaries. This increased independence between sediment concentration and discharge may be attributed to the heterogeneity in these relatively large catchments.

This emergent stationarity explains the linear correlation between area-specific sediment yield and runoff depth reported in a small sub-watershed in a hilly area of the Loess Plateau (Zheng et al., 2013). Considering the sediment concentration to be constant, the variation in yield is solely dominated by streamflow, resulting in the observed linear discharge–yield relationship. Similar stationarity in sediment concentration has also been found in arid watersheds in Arizona, US (Gao et al., 2013), where the sediment concentration becomes homogeneous among watersheds when their drainage area is larger than 0.01 km². The difference in threshold for the emergence of approximately discharge-invariant concentration between the YRB and watersheds in Arizona, US, is probably due to the differences in catchment characteristics, i.e., vegetation type and coverage, terrestrial structure, soil properties, etc.

Our analysis shows that mean annual discharge (Q_m) is a better indicator of the correlation between water and sediment transport than drainage area, although the last parameter has been used traditionally. Despite the heterogeneity, both the coupling between Q–C and the concentration sensitivity to variation in streamflow decreases with Q_m. A closer inspection reveals useful insights. At gauges with smaller values of Q_m, discharge is the dominant factor in sediment transport: an increment in annual discharge is amplified in the increment of sediment concentration (α_QC>1) (i.e., Gauge 808, 812 in Fig. S4). However, as Q_m increases, variation in streamflow is more weakly reflected in variation in sediment concentration, even though annual mean concentration still correlates with annual discharge, (i.e., Gauge 806 in Fig. S4). As Q_m continues to increase, sediment concentration becomes almost invariant to discharge, suggesting that the dominant factor of sediment transport has shifted from the discharge to something else.

4 The vegetation impact: a bell shape

To further explore the potential cause of this emergent stationarity, we analyzed the vegetation data (NDVI) from 44 of the gauges located on eight major tributaries of the YR (Fig. S1). Our analysis shows that this declining sensitivity in concentration at an annual scale (α_QC) is negatively related to vegetation impact across the gauges (Fig. 2).

Figure 2 Scatter plots between the maximum NDVI and slope in the Q–C regression at an annual scale (α_QC) from the 44 study gauges.

For gauges with limited vegetation establishment in their drainage area, the variation in discharge is amplified in sediment transport (α_QC>1). The larger the discharge is in a specific year, the more sediment is eroded and mobilized per cubic meter. This dominance of discharge is weakened when vegetation density and coverage increase. Despite the larger sediment carrying capacity of larger discharge, sediment concentration is reduced, probably due to the protection vegetation offers against erosion. As maximum NDVI increases, sediment concentration becomes less and less coupled with discharge at an annual scale. When the vegetation density is sufficiently high, sediment concentration is nearly stable in spite of the variation in discharge, since the dense vegetation coverage protects soil from erosion and traps sediment. That is, the emergent stationarity in sediment concentration corresponding to the variation in discharge at gauges with large Q_m can be attributed to the dampened dominance of discharge due to the increasing impact of vegetation retardation.

Figure 3 Scatter plot of annual mean concentration and maximum NDVI: (a) at 44 study gauges between 1982 and 1986, where the dots are color-coded by the slope in the Q–C regression (α_QC) at each gauge; and (b) at seven gauges with both data from the years 1982 to 1986 (blue dots) and the years 2008 to 2012 (red dots). The R² for the two fits is 0.6 and 0.44 respectively with a p value <0.001 for both of them.

To further confirm the vegetation impact on sediment transport, we derived the plot between maximum NDVI and mean concentration at an annual scale in Fig. 3a. As we can see, the annual mean sediment concentration follows a bell-shaped correlation with vegetation establishment, with a peak concentration at a value of maximum NDVI of around 0.36. On the falling limb of this bell curve, as NDVI increases, both sediment concentration and α_QC decrease consistently. That is, both the value of concentration and its sensitivity to streamflow variation decline with increasing vegetation index on the falling limb. To confirm this impact of vegetation resistance, we also examined the relationship between sediment concentration and other catchment characteristic like dominant soil type. No significant correlation was observed. Although there could still be other factors not considered here that contributed to the decline in sediment concentration, it is undoubted that vegetation is one of the most influential factors of sediment reduction and can be used as a good indicator of the soil erosion and sediment transport in the YRB.

On the rising limb, however, both the value of concentration and its sensitivity to streamflow variation increase with increasing vegetation index. Most gauges have values α_QC larger than 1, except one gauge with an extremely small maximum value of NDVI. For these gauges, on the rising limb, vegetal cover is still low in an absolute sense despite increasing NDVI. Sediment concentration is mainly dominated by discharge: fluctuations in streamflow are amplified in concentration (α_QC>1). The only gauge with a value of α_QC smaller than 1 is gauge Hanjiamao (HJM) at the Wuding River. Although the annual precipitation and discharge at HJM is similar to other gauges along the Wuding River, the annual mean sediment concentration is much smaller. This is because of the extremely high baseflow contribution in discharge at HJM, which is around 90 %, thanks to very intensive check-dam construction there (Dong and Chang, 2014). Since sediment in the YRB is mostly transported during large flow events during the summer, smaller flow events are not capable of transporting significant sediment loads at HJM.

Figure 4 Illustration of the correlation between vegetation and sediment erosion, retardation, and the resulting sediment concentration in the YRB. Since vegetation usually increases with discharge, with the rise in discharge, sediment eroded and delivered by streamflow increases rapidly, while the retardation from vegetation is limited at the beginning and increases fast afterwards. This non-synchronous impact on sediment transport leads to the bell shape correlation between sediment concentration and vegetation.

In general, we can conclude that sediment transport is mainly dominated by discharge when the vegetation index is low. With increasing NDVI, the impact of vegetation grows slowly at first and accelerates after the maximum NDVI exceeds 0.36. Eventually, the effect of NDVI takes over the dominance of streamflow and attenuates the variation in sediment concentration (Fig. 4). The nonlinear impact of vegetation in regard to resistance of sediment to erosion is consistent with previous findings (Rogers and Schumm, 1991; Collins et al., 2004; Temmerman et al., 2005; Corenblit et al., 2009). When the vegetation index level is low, its resistance to soil erosion develops slowly as vegetation grows and expands (Rogers and Schumm, 1991), and the capability of vegetation to trap sediment is reduced when submerged by flood (Temmerman et al., 2005) or overland flow. Therefore, for catchments with limited vegetation establishment, the coverage of vegetation is not sufficient to trap sediment, nor is the vegetation able to protrude from the water level during the extreme flow events that transport most of the sediment. Sediment transport in these catchments is usually dominated by discharge. As NDVI increases, vegetation becomes much more capable as an agent of erosion protection and sediment settling (Jordanova and James, 2003; Corenblit et al., 2009). With the compensation from vegetation retardation, sediment and discharge become more and more decoupled as discharge increases, so concentration is nearly invariant to increasing discharge. The transition point in maximum NDVI (around 0.36) is where the increment in vegetation reduction balances with the incremental increase in water erosion. When the capability of vegetation retardation catches up with streamflow erosion, the net soil loss becomes negligible, a condition commonly observed in well-vegetated regions.

5 Validation of the bell shape across time and space

In 1999, a large-scale ecosystem restoration project, the Grain-for-Green project, was launched in the YRB for soil conservation (Lv et al., 2012). It has substantially improved vegetation coverage after a decade of implementation (Sun et al., 2015). To validate our hypothesis gain from the early 1980s, we applied a similar analysis to the annual flow and sediment data as well as daily NDVI data at seven gauges located at the outlets of major tributaries from 2008 to 2012 (Fig. S1 green stars). This is the period subsequent to the initiation of the Grain-for-Green project. We have excluded the years right after the implementation of the Grain-for-Green project, when there was an initial drastic change in vegetation coverage and sediment erosion and transport processes.

As we can see from Fig. 3b, there is a significant increase in maximum NDVI for all seven catchments and considerable reduction in mean sediment concentration. This improvement is consistent with the previous report that the Grain-for-Green project has made a remarkable achievement in regard to soil conservation in the YRB (Chen et al., 2015). Comparison of the relationship between sediment concentration and maximum NDVI in the early 1980s and around 2010 shows that the bell shape relationship is sustained even after drastic and significant anthropogenic alteration of the land use and land cover across the whole YRB. Although the vegetation coverage has improved significantly at all seven comparison gauges due to the ecosystem restoration policy, and therefore effectively moderated sediment erosion, the bell shape relationship between maximum NDVI and mean concentration is sustained.

Similar bell shape relationship was also found for the multi-year mean annual precipitation and sediment yield observed in the United States (Langbein and Schumm, 1958). The data used in the analysis of Langbein and Schumm (1958) were collected in the 1950s from more humid and vegetated catchments with limited human intervention, opposite to the YRB where the climate is arid and semi-arid, vegetation coverage is low, and human activities are intensive. However, a similar bell shape was still observed between sediment yield and precipitation. Given the limited anthropogenic activities in these catchments, vegetation growth is probably correlated with annual precipitation due to its adaption to climate, as in other US catchments (Fig. S6). Thus it is likely that a bell shape correlation between vegetation and sediment yield would be found at these US catchments as well. This suggests that the bell shape correlation between vegetation and sediment concentration is not only observed in the YRB with intensive human intervention, but could also be valid outside it. More analyses are needed to test this relationship in other catchments outside the YRB for its universality.

6 Implications and conclusion

Our analysis shows that across the YRB both the correlation between Q and C and the magnitude of sediment response to the variation in streamflow decrease with Q_m. When Q_m is sufficiently large (i.e., >60 mm yr⁻¹), sediment concentration reaches a stationary (constant) state at an annual scale. The emergent stationarity at gauges with large Q_m is related to the shift of dominance from discharge to vegetation. Because of the slow development of vegetation resistance with increasing discharge for small discharges, discharge dominates the soil erosion and sediment transport process until the maximum NDVI exceeds a threshold (0.36 for this study) at which the parameter governing concentration transits from streamflow erosion to vegetation retardation.

Our findings of the emergent stationarity in sediment concentration and the shift of the dominant mechanism governing the Q–C relation have important implications for water and sediment management at the watershed scale. Our study indicates that, for the gauges with relatively large discharge, the annual mean concentration can be approximated as a constant over a large range of discharges. Thus the estimation of sediment yield can be simply inferred from a simulation of streamflow. First-order estimates of sediment yield for scientific or engineering purposes can be obtained by multiplying the estimated discharge by a constant sediment concentration estimated based upon the vegetation index. The correlation between vegetation and sediment concentration will also be useful for the design of the ongoing ecosystem restoration program known as the Grain-for-Green project. The bell-shaped correlation between maximum NDVI and sediment concentration provides a quantitative way to estimate the potential change in sediment concentration associated with proposed ecosystem restoration planning schemes at and near each tributary. This can help guide land use management so as to allocate the sediment contribution from each of the upstream tributaries in a way that maintains the balance between erosion and deposition in the lower YR.

It is important to collect more data from the current decade (i.e., after the substantial ecosystem restoration) to further validate our findings in regard to emergent stationarity and vegetation impact at more gauges in the YRB. It will be helpful if we could examine our findings in other watersheds worldwide with different climate and vegetation types. Although humid regions are usually considered to be well vegetated, studies show that there could still be erosion issues in these areas due to topographic gradient, precipitation intensity, and soil properties. (Holz et al., 2015). Analysis with more field measurements could also help explain the threshold discharge of the emergent stationarity. Numerical simulations as well as long-term measurements on the soil properties are also needed to further explain the physical mechanism of vegetation retardation: how it develops its impact on soil erosion and sediment transport by changing soil properties and other topographic characteristics during its growth and spread.