Estimated Grass Grazing Removal Rate in a Semiarid Eurasian Steppe Watershed as Influenced by Climate

Grazing removal rate of grasses needs to be determined for various climate conditions to address eco-environmental concerns (e.g., desertification) related to steppe grassland degradation. The conventional approach, which requires survey data on animal species and heads as well as grass consumption per individual animal, is too costly and time-consuming to be applied at a watershed scale. The objective of this study was to present a new approach that can be used to estimate grazing removal rate with no requirement of animal-related data. The application of this new approach was demonstrated in a Eurasian semiarid typical-steppe watershed for an analysis period of 2000 to 2010. The results indicate that the removal rate tended to become larger, but its temporal variation tended to become smaller, from the upstream to downstream. Averaged across the watershed, the removal rate ranged from 63.9 to 401.0 g DM m ́2 (or 22.4 to 60.9%) during the analysis period. As expected, the removal rate in an atmospherically wetter year was higher than that in an atmospherically drier year. Nevertheless, none of the eleven analysis years had a removal rate higher than the threshold value of 65%, above which the risk of grassland degradation would become much greater.


Introduction
Steppe grasslands occupy about 8% of the earth's terrestrial surface and are now considered the most altered and beleaguered ecosystem on the planet [1].Over 40% of the global grasslands have been somewhat altered from their indigenous state, including more than 70% of the Inner Mongolia steppe grasslands of China [2], which are part of the largest and most characteristic Eurasian grassland (or Great Steppe).Here, the degradation has caused serious environmental and ecological problems such as more frequent and damaging hydrologic extremes (i.e., floods and droughts), desertification, dust storms, and even commodity scarcity.These problems in turn will likely threaten the sustainability of "grassland agriculture," a system of agriculture in which major emphasis is placed on grasses, legumes, and other fodder or soil-building crops [3].The main factors causing degradation are overgrazing, cultivation, overdevelopment, and climate change [4][5][6][7][8], which have altered the natural hydrology [9] and led to the erosion of the 10 to 20 cm calcic castanozem topsoil.This topsoil is vital for efforts to sustain grasslands because it is loose and has a plentiful supply of humus/organic matter.Policy choices to reduce or reverse grassland degradation are often made with only a vague understanding of causative complexity among steppe hydrology, topsoil erosion, and grassland degradation, however, which can limit efforts to manage, protect and/or restore steppe grasslands.
In order to fill this knowledge gap, the grazing removal rate of grasses needs to be determined for various climate conditions as characterized by precipitation and evapotranspiration.The conventional approach consists of that: (1) an extensive survey is conducted household by household to obtain data on species and heads of all grass-grazing animals (i.e., cows, horses, and sheep); (2) for a given animal species, the grazing rate was computed as the multiplication of the heads of this species and the dry mass (DM) of grass eaten by one such animal per year; and (3) the total grass removal rate in an area of interest is computed as the summation of the grazing rates of all animal species.Such a conventional approach can become too costly and time-consuming to be implemented for multiple years at a watershed scale as required to effectively manage steppe grasslands in practice.This may be an important reason why few watershed-scale studies on grazing removal rate have been reported in existing literature.The objective of this study was to present an innovative method that can be used to accurately estimate grazing removal rate at a watershed scale using remote sensing images and a total net primary production (NPP) prediction model.Hereinafter, NPP is the summation of aboveground (i.e., stems and leaves) and belowground (i.e., roots) productions.One advantage of this method over the conventional approach is that this method does not require data on animals.The rationale of this method is that NPP can well reflect the physical link between biosphere and climate system through the global cycling of carbon, water and nutrients, and can be reduced by grazing as a primary mechanic mechanism.This study applied and tested this method in a Eurasian semiarid typical-steppe watershed to be described in Section 2.1.

The Study Watershed
The 5350 km 2 Balagaer River watershed (centroid: 117 ˝36 1 E, 44 ˝35 1 N), located in the northeast Inner Mongolia Autonomous Region of China (Figure 1), was selected for this study.This watershed is almost uniformly covered by typical steppe grasses (Figure 2), with the predominant species of Leymus chinensis and Stipa grandis.The grasses, which are solely fed by rainwater and dependent on naturally-available nutrients (e.g., animal manures), have a root depth of up to 60 cm [10].Rarely harvested, most of the grasses are mainly removed by semi-nomadic-style (i.e., half-year) grazing.The elevation of the watershed varies from 980 to 1876 m above mean sea level, with a mean topographic gradient of 0.09.The soils are classified by the Food and Agriculture Organization of the United Nations [11] into four classes, namely I-Mo-2c, Kh1-2b, I-K-2c, and Zm2-2/3a (Figure 1), each of which is further subdivided into two (i.e., upper and lower) layers.The first three types of soils are more permeable than the fourth.The upper layer is from ground surface to the 30 cm depth, while the lower layer is from the 30 to 100 cm depth.Across the watershed, the saturated hydraulic conductivity (K s ) varies from 1.5 to 107.0 mm¨h ´1, with a mean of 28.5 mm¨h ´1.The loamy soils are composed of less than 30% clay particles (diameter less than 0.002 mm) and more than 70% sand (diameter from 0.05 to 2.0 mm) and silt (diameter from 0.002 to 0.05 mm) particles [12].
The watershed receives 170 to 615 mm precipitation annually, with a mean of 335 mm.Most of the precipitation falls between July and September as rain, and between October and January as snow.However, the watershed has an annual potential evapotranspiration (E 0 ) of 1165 mm or higher, which is much larger than the mean annual precipitation, so the watershed has an arid and semiarid climate (i.e., is water-limited).The annual mean daily average air temperature is 1.2 ˝C, with a maximum daily temperature of up to 37.5 ˝C in summer and a minimum daily temperature of as low as ´38.5 ˝C in winter.The watershed has an annual mean daily average wind speed of about 15 km¨h ´1, with between 28 and 148 windy days in any given year.On a windy day, the maximum wind speed can reach 125 km¨h ´1.The watershed has an annual mean discharge of 0.33 m 3 ¨s´1 , and the major water users are agriculture, animal husbandry, domestic, and power industry [13].

Data and Preprocessing
Daily data on precipitation, air temperature, relative humidity, wind speed, and solar radiation at West Ujimqin Climate Station 54012 (Figure 1) were downloaded from the National Meteorological Information Center website [14] for a record period of 1 January 2000 to 31 December 2010.This study presumed that precipitation was spatially uniform across the watershed for three reasons.First, there was no any other climate station within or adjacent to the watershed, where long-term observations were available.Second, given that West Ujimqin Climate Station is located near the centroid of the watershed, the station could likely represent the primary climate conditions for the entire drainage area.Thirdly, using the data on precipitation measured from 2010 to 2012 at 15 raingauges, which were installed and maintained by the authors across the watershed, Luo et al. [15] found that the precipitation exhibited a minimal spatial variation.In fact, the sufficient spatial coverage of climate stations and/or raingauges has been, and will continue to be, a challenge for any hydrology-related studies provided that the insufficient coverage tends to under-represent the possible spatial heterogeneity of precipitation.Global efforts are underway to use modern technologies (e.g., next generation radar) to acquire precipitation data at a kilometer spatial

Data and Preprocessing
Daily data on precipitation, air temperature, relative humidity, wind speed, and solar radiation at West Ujimqin Climate Station 54012 (Figure 1) were downloaded from the National Meteorological Information Center website [14] for a record period of 1 January 2000 to 31 December 2010.This study presumed that precipitation was spatially uniform across the watershed for three reasons.First, there was no any other climate station within or adjacent to the watershed, where long-term observations were available.Second, given that West Ujimqin Climate Station is located near the centroid of the watershed, the station could likely represent the primary climate conditions for the entire drainage area.Thirdly, using the data on precipitation measured from 2010 to 2012 at 15 raingauges, which were installed and maintained by the authors across the watershed, Luo et al. [15] found that the precipitation exhibited a minimal spatial variation.In fact, the sufficient spatial coverage of climate stations and/or raingauges has been, and will continue to be, a challenge for any hydrology-related studies provided that the insufficient coverage tends to under-represent the possible spatial heterogeneity of precipitation.Global efforts are underway to use modern technologies (e.g., next generation radar) to acquire precipitation data at a kilometer spatial

Data and Preprocessing
Daily data on precipitation, air temperature, relative humidity, wind speed, and solar radiation at West Ujimqin Climate Station 54012 (Figure 1) were downloaded from the National Meteorological Information Center website [14] for a record period of 1 January 2000 to 31 December 2010.This study presumed that precipitation was spatially uniform across the watershed for three reasons.First, there was no any other climate station within or adjacent to the watershed, where long-term observations were available.Second, given that West Ujimqin Climate Station is located near the centroid of the watershed, the station could likely represent the primary climate conditions for the entire drainage area.Thirdly, using the data on precipitation measured from 2010 to 2012 at 15 raingauges, which were installed and maintained by the authors across the watershed, Luo et al. [15] found that the precipitation exhibited a minimal spatial variation.In fact, the sufficient spatial coverage of climate stations and/or raingauges has been, and will continue to be, a challenge for any hydrology-related studies provided that the insufficient coverage tends to under-represent the possible spatial heterogeneity of precipitation.Global efforts are underway to use modern technologies (e.g., next generation radar) to acquire precipitation data at a kilometer spatial resolution [16].However, such high-resolution data are just available for a few experiment sites and can have large uncertainties [17].
In addition, a soil classification GIS (geographic information system) map and its associated attribute table were downloaded from the FAO website [18].The GIS map subdivides the watershed into 22 map units (not shown), each of which is assigned one of the four classes mentioned in Section 2.1, as shown in Figure 1.For each soil map unit and layer, the attribute table defines its soil texture (i.e., percent sand, silt, and clay particles) and organic matter content.In this study, the information presented by the attribute table was used to estimate the intrinsic soil-water properties (e.g., permanent wilting point ψ, field capacity θ fc , and saturated soil moisture θ sat ) [19].ψ is the minimal point of soil moisture the plant requires not to wilt, θ sat , on the other hand, is the soil moisture when all pores in the soil are filled with water.θ fc is the soil moisture after the soil is saturated and then allowed to drain by water gravity until the drainage ceases.Further, a 30-m digital elevation model (DEM), which is a 3D representation of the watershed's terrain surface with a spatial resolution of 30 m, was downloaded from the International Scientific Data Service Platform website http://datamirror.csdb.cnand used to delineate the boundary of the study watershed and its subbasins.As a result of the delineation, the watershed was subdivided into 27 subbasins (Figure 1).
Further, data on annual NPP, in g C m ´2, were extracted using the ENVI ® 4.7 software package from the 1-km MODIS (Moderate Resolution Imaging Spectroradiation) 17A3 images for the period of 2000 to 2010.The images, downloaded from the U.S. Geological Survey (USGS) Land Processes Distributed Active Center (LPDAC) website [20], present values of annual NPP [21][22][23][24].Firstly, in ArcGIS ® 10, the images that cover the study watershed were merged into a single image, which was then projected to the coordinate system of the Universal Transverse Mercator (UTM) zone 50 N.Secondly, the projected image was clipped to the spatial extent of the watershed and then overlaid with the subbasin map delineated above.Thirdly, for each year, the NPP values for the 1-km cells that are included within a subbasin were arithmetically averaged to get the NPP value of this subbasin, and then the NPP values for the 27 subbasins were area-weighted averaged to get the NPP value of the study watershed as a whole.As a result, 28 NPP time series (one for the watershed and 27 for the subbasins) were generated.Finally, for a given subbasin, the NPP values for the analysis years (2000 to 2011) were arithmetically averaged to get the annual mean NPP of this subbasin.

Description of the NPP-Prediction Model
This study modified and used an NPP-prediction model developed by [13].This empirical model estimates NPP, in g DM m ´2, as: NPP " 100 ¨RDI ¨E ¨exp " ´p9.87 `6.25 ¨RDIq 0.5 ı where RDI (-) is the dimensionless climatic dryness index; E (mm) is the annual actual evapotranspiration; and 100 is a unit conversion factor.RDI is computed as [25]: where E 0 (mm) is the annual potential evapotranspiration; and P (mm) is the annual precipitation.E is estimated using a modified theoretical Fu's solution [26][27][28] to the hypothesis of [29].The Budyko hypothesis describes the annual water balance as a function of available water and energy, and it has been tested all over the world and widely used to estimate E [30][31][32].However, because the original Fu's solution (used in [13]) assumes that E ď P, which can be invalid when P is small in a water-limited environment [33], which is the case of the Balagaer River watershed, this study developed a modified solution by introducing a soil-moisture enhancement factor to take into account the evapotranspiration portion beyond P. Albeit, both the modified and original Fu's solution share a same assumption of negligible water storage change within the analysis time interval.This assumption is usually valid for arid and semiarid areas at a large (e.g., annual) time interval [34].The modified Fu's solution is expressed as: where α sm (ě1) (-) is the soil-moisture enhancement factor, representing the compensation of soil moisture to the portion of E beyond P; ω (-) has a range of [1, +8) and is a parameter that represents combined effects of subbasin characteristics (namely vegetation physiology, soil properties, and topography) on water balance.Yang et al. [30] presented a formula for ω.This formula (used in [13] considers soil properties and topography, but it does not have any variable representing vegetation physiology.Based on basic principles of hydrology [35] and findings of [36], a vegetation canopy can enhance evapotranspiration, that is, E tends to become larger for an area with a vegetation canopy than for the same area without the canopy.The enhancement depends on vegetation physiologic characteristics, namely species, stem height, leaf structure, and root structure.With this regard, this study introduced a vegetation-specific enhancement factor as a multiplier into the formula of [30], resulting in a new formula for ω expressed as: where D w (-) is the number of wet days with nonzero precipitation; K s (mm¨day ´1) is the saturated hydraulic conductivity; S max = (θ fc ´ψ)¨d r (mm) is the available water capacity of root zone; d r (mm) is the root depth; β ( ˝) is the landscape slope angle; and α veg (ě1) (-) is the vegetation-specific enhancement factor.Because the Penman-Monteith formula [37] is physically based and the best predictor of potential evapotranspiration (PET) among the existing evaporation models [38,39], it was chosen to estimate daily PET (mm¨day ´1), which was summed spanning a year to get E 0 in this NPP-prediction model.The formula can be expressed as: PET " 0.408 ¨∆ ¨pR n ´Gq `γ ¨900 T air `273 ¨u2 ¨pe s ´ea q ∆ `γ ¨p1 `0.34 ¨u2 q (5) where ∆ (kPa¨˝C ´1) is the slope of saturation vapor pressure curve; R n (MJ m ´2¨day ´1) is the net solar radiation at vegetation surface; G (MJ m ´2¨day ´1) is the soil heat flux density; γ (kPa¨˝C ´1) is the psychrometric constant; T air ( ˝C) is the mean air temperature at 2 m height above ground surface; u 2 (m¨s ´1) is the wind speed at 2 m height above ground surface; e s (kPa) is the saturation vapor pressure; and e a (kPa) is the actual vapor pressure.These parameters can be determined using the procedure recommended by [11].T air is computed as: where T min ( ˝C) and T max ( ˝C), respectively, are the daily minimum and maximum air temperatures.
∆ is computed as: ∆ " p0.00815 ¨Tair `0.8912qWhen wind speed is measured at a height different from 2 m above the ground, it can be converted to u 2 using a logarithmic profile formula found in [41,42].R n can be directly measured or estimated based on the extraterrestrial solar radiation and air temperatures [35,43].

Parameterization of the NPP-Prediction Model
For each soil map unit, the soil-water parameters, includingψ, θ fc , θ sat , and K s , were estimated using the soil-plant-air-water (SPAW) model [18].The model uses a set of empirical equations to formulate the soil-retention curves in terms of the soil texture.The curves include the pressure-moisture curve, which describes the continuous relationship between capillary pressure and soil-water moisture content, and the conductivity-moisture curve, which describes the relationship between hydraulic conductivity and soil-water moisture content.In general, a pressure-moisture curve consists of three segments [44,45]: the first segment is continuous and nonlinear for capillary pressures between 1500 and 10 kPa, the second segment is linear for capillary pressures between 10 kPa to the air-entry capillary pressure, and the third segment has constant water content for capillary pressures smaller than the air-entry capillary pressure.In contrast, the conductivity-moisture curve is continuous and nonlinear from the saturation moisture content (θ sat ) to near air dry.Those empirical equations, described in detail in [46], have been proved to be valid for a wide range of textures.The inputs of the SPAW model include percentages of sand and clay particles, and the average organic matter content within a soil map unit of interest.In addition, the ArcGIS ® Hydrology extension was used to process the DEM to delineate the boundary of the Balagaer River watershed and its subbasins.As a result of the delineation, the watershed was subdivided into 27 subbasins (Figure 1), each of which was assumed to have a uniform overland gradient (i.e., single value for β).The extension automatically calculated β for the subbasins from the elevation values presented by the DEM.
Further, the subbasin map and the soil map were overlain in ArcGIS ® to determine the areal proportion of a given soil map unit that is included in each of the subbasins.Afterward, for each of the four parameters (i.e., ψ, θ fc , θ sat , and K s ), the area-weighted average of its values responding to the map units that are included in a subbasin was computed and taken as the value of this parameter for this subbasin.Moreover, this study assumed that d r (=60 cm), P atm (=101.3kPa), α sm , and α veg are spatiotemporally invariant.α sm and α veg were empirically adjusted until E P for the watershed and its subbasins falls between 0.85 to 1.35 for all eleven analysis years.This range of E P was determined by [33] based on field measurements and the FAO-56 (Food and Agriculture Organization Irrigation and Drainage Paper No. 56) dual crop coefficient (DCC) method [47,48] in the Xilin River Basin (115 ˝32 1 to 117 ˝36 1 E, 43 ˝26 1 to 44 ˝39 1 N), which is adjacent to the Balagaer River watershed.

Analysis of Grazing Removal Rate
For a given year, the value predicted by the model (Equation ( 1)) is the NPP with no influence of grazing and was thus presumed as the "natural" NPP, whereas the value presented by the MODIS data reflects the ground truth and was considered as the "actual" NPP.With a typical carbon content of 70% in dry leaves/stems plus roots of steppe grasses (including Leymus chinensis and Stipa grandis) [49], the MODIS-presented actual NPP (in g C m ´2) was divided by this carbon content to have consistent units with the mode-predicted natural NPP (in g DM m ´2).For the simplicity, by assuming that grazing and climate independently (but not interactively) influenced grass production, this study computed grass removal rate as difference of natural NPP less actual NPP.Although previous field experiments (e.g., [50][51][52]) found that the grass production in a dry year was mainly controlled by precipitation and temperature while the grass production in a wet year was primarily affected by grazing, [53] postulated that grazing and climate can somewhat interactively influence grass production and that such an interactive influence can either be positive or negative, depending on a number of not-easy-to-be-quantified factors (e.g., grazing intensity and pattern).Thus, our simple method may overestimate the grass removal rate for one year, while it can underestimate the grass removal rate for another year.However, for a long-term average, the overestimations and underestimations can be expected to be canceled out, making the method be a good predictor of annual mean NPP.
The computation, which was done on a yearly basis for each of the 27 subbasins as well as for the study watershed as a whole, resulted in 28 time series of removal rate in g DM m ´2, which in turn were divided by the corresponding 28 time series of natural NPP in g DM m ´2 to derive another 28 time series of removal rate in percent.Totally, 112 (4 ˆ28) time series, each of which has a length of N = 11 years, were obtained: 28 for natural NPP, 28 for actual NPP, 28 for removal rate in g DM m ´2, and 28 for removal rate in percent.Please note that four of the 112 time series were for the watershed and the remaining 108 time series were for the subbasins.In addition, 58 time series (with a same length of N = 11 years) of climate factors were also derived: one for P, 28 for E, 28 for E P , and one for E 0 P .Further, across the 27 subbasins, the maximums and minimums of natural NPP, actual NPP, removal rate in g DM m ´2, and removal rate in percent, were computed, resulting in eight more time series of extreme values, each of which has a length of N = 11 years.
Firstly, the eight time series of extreme values and the four time series for the watershed were plotted to visually examine the spatiotemporal variations of NPP and removal rate.Secondly, at the watershed scale, the removal rate in percent was plotted versus each of the four climate factors (i.e., P, E, E P , and E 0 P ) to visually examine overall influences of climate on grazing removal rate.Finally, at both the watershed and subbasin scales, Pearson correlation coefficients [54] were computed between removal rate (either in g DM m ´2 or percent) and each of the four climate factors, and then compared with a critical value r c (Equation ( 12)) for a significance level of α = 0.05 to identify which climatic condition (e.g., drier versus wetter) might be more closely associated with a higher grazing removal rate.A climate factor that gives a Pearson correlation coefficient (in absolute) of greater than r c was judged to have a significant (at α = 0.05) influence on grass removal.
where N (years) is the length of time series; and t α (-) is the critical value of a t-distribution with N ´2 freedom and at a two-tail exceedance probability of α.Here, r c " " 0.60.

The Parameterized Model
The sizes of the delineated subbasins ranged from 3 to 458 km 2 , with a mean of 198 km 2 .As expected, the subbasins that are directly drained by the lower reach of the Balagaer River Water 2016, 8, 339 8 of 18 (e.g., subbasins 3 and 6) were determined to have a smaller landscape slope angle than the subbasins near the watershed upper boundary (e.g., subbasins 20 and 26) (Figure 1 and Table 1).The calculated subbasin slope angles vary from 1.29 ˝to 8.67 ˝, with an overall watershed slope angle of 4.17 In addition, subbasin 20 was determined to have largest values for ψ, θ fc , and θ sat , whereas subbasin 9 was determined to have smallest values for ψ and θ fc .Subbasins 7 and 8 were determined to have a smallest value for θ sat .At the watershed level, the soils were determined to have an overall retention storage of θ fc ´ψ = 0.138 and an overall detention storage of θ sat ´θfc = 0.183.Further, the subbasins that are directly drained by the middle and upper reaches of the Balagaer River (e.g., subbasins 9 and 12) were determined to have higher values for K s (i.e., be more permeable) than the subbasins that are directly drained by the lower reach.Moreover, for the watershed, the average saturated hydraulic conductivity was determined to be 22.69 mm¨h ´1.3)).The parameters of β, ψ, θ fc , θ sat , and K s are used in Equation ( 4); ‡ See Figure 1.
The vegetation-specific enhancement factor was taken as α veg = 2.5, which is equal to the long-term average leaf area index (LAI) of matured steppe grasses [42], while the soil-moisture enhancement factor was adjusted to be α sm = 1.15, which gave E P " 0.9 to 1.3 (Figure 3), depending on the subbasin and/or analysis year of interest.Overall, while E tended to increase with P by following a quadratic function (coefficient of determination R 2 = 0.73), it varied greatly from one subbasin to another for a given P. In contrast, E P was found to be barely dependent on P, as indicated by the large scattering of the points and a much smaller R 2 = 0.18.
Water 2016, 8, 339 9 of 18 P on the subbasin and/or analysis year of interest.Overall, while E tended to increase with P by following a quadratic function (coefficient of determination R 2 = 0.73), it varied greatly from one subbasin to another for a given P. In contrast, P E was found to be barely dependent on P, as indicated by the large scattering of the points and a much smaller R 2 = 0.18.

The Estimated Removal Rate of Grasses
For the study watershed as a whole, the NPP could have varied from 284.7 g DM m −2 in 2007 to 687.4 g DM m −2 in 2008 if there were no grazing, whereas the actual NPP was determined to be 220.8g DM m −2 in 2007 to 324.7 g DM m −2 in 2003 (Figure 4a).This indicated that 63.9 to 401.0 g DM m −2 grasses were grazed (Figure 4b), which is equivalent to 22.4% to 60.9% of the non-grazing grass production (i.e., the predicted natural NPP) (Figure 4c).
Across the study watershed, for a given analysis year, the natural NPP was predicted to be 282.6 g DM m −2 in one subbasin but as high as 704.0 g DM m −2 in another subbasin, while the actual NPP was determined to be 161.6 g DM m −2 in one subbasin but 371.6 g DM m −2 in another subbasin (Figure 4a).Grazing removed 8.7 to 437.9 g DM m −2 at the subbasin scales (Figure 5a), which is equivalent to 3.1% to 74.1% of the non-grazing grass production (Figure 5b).Overall, the subbasins downstream of the West Ujimqin banner (near the Baiyingwula flow station shown in Figure 1) had a significantly (p-value = 0.01 < α = 0.05) higher grazing rate than the subbasins upstream of the banner, as indicated by a paired student-t (unequal variances) test [54].The null hypothesis of the test is that the mean grazing rate of the upper portion of the watershed was not different from the mean grazing rate of the lower portion.However, the grazing rates in the downstream subbasins had a larger variation from year to year than those in the upstream subbasins.The subbasins with maximums and minimums for the analysis years are presented in Table 2.

The Estimated Removal Rate of Grasses
For the study watershed as a whole, the NPP could have varied from 284.7 g DM m ´2 in 2007 to 687.4 g DM m ´2 in 2008 if there were no grazing, whereas the actual NPP was determined to be 220.8g DM m ´2 in 2007 to 324.7 g DM m ´2 in 2003 (Figure 4a).This indicated that 63.9 to 401.0 g DM m ´2 grasses were grazed (Figure 4b), which is equivalent to 22.4% to 60.9% of the non-grazing grass production (i.e., the predicted natural NPP) (Figure 4c).
Across the study watershed, for a given analysis year, the natural NPP was predicted to be 282.6 g DM m ´2 in one subbasin but as high as 704.0 g DM m ´2 in another subbasin, while the actual NPP was determined to be 161.6 g DM m ´2 in one subbasin but 371.6 g DM m ´2 in another subbasin (Figure 4a).Grazing removed 8.7 to 437.9 g DM m ´2 at the subbasin scales (Figure 5a), which is equivalent to 3.1% to 74.1% of the non-grazing grass production (Figure 5b).Overall, the subbasins downstream of the West Ujimqin banner (near the Baiyingwula flow station shown in Figure 1) had a significantly (p-value = 0.01 < α = 0.05) higher grazing rate than the subbasins upstream of the banner, as indicated by a paired student-t (unequal variances) test [54].The null hypothesis of the test is that the mean grazing rate of the upper portion of the watershed was not different from the mean grazing rate of the lower portion.However, the grazing rates in the downstream subbasins had a larger variation from year to year than those in the upstream subbasins.The subbasins with maximums and minimums for the analysis years are presented in Table 2.

The Influences of Climate on Grazing Removal Rate of Grasses
The grass removal rate tended to significantly increase with P, E, or E P , whereas it tended to marginally decrease with E 0 P (i.e., atmospheric dryness) (Table 3 and Figure 6).A year with a greater E 0 P is drier, which can restrict vegetation growth resulting in less available grasses for grazing.In contrast, a year with a smaller E 0 P tends to be wetter, which can increase vegetation growth leading to more available grasses for grazing.Among these four climate factors, E was found to be most positively related to grass production.This is because E reflects the overall natural conditions for vegetation photosynthesis, with a higher E indicating a better condition; and vice versa.In terms of removal rate in percent, the downstream subbasins were less likely influenced by P and E 0 P than the upstream subbasins, as indicated by that more downstream subbasins have a statistically insignificant Pearson correlation coefficient (Table 3).However, such a contrast between the upstream and downstream subbasins was not found in terms of removal rate in g DM m ´2.This indicates that the downstream subbasins might have a consistently high removal percentage regardless of natural grass production and year (i.e., had a proportionally increased grazing intensity with natural grass production).Herein, it is assumed that grazing intensity and removal rate are positively correlated: a higher grazing intensity corresponds to a larger removal rate, and vice versa.The grazing intensities in the upstream subbasins might have large variations on a yearly basis.

Discussion
The method (Equation (1) through Equation ( 11)) can be parameterized using easily-available data on topography, soil, land use and land cover, and climate.For an area of interest, the method predicts the natural NPP, whereas the MODIS images present the actual NPP.The removal rate can be estimated as the difference of the natural NPP less the actual NPP.In comparison to the conventional survey-based approach, the method is universally applicable and does not require animal-related data (e.g., species and heads).The method can be applied to areas with a small (or large) size and an unknown simple (or complex) grazing pattern.Leriche et al. [53] proposed a model for evaluating the short-term effect of grazing on grass NPP.This model requires inputs of grazing intensity and pattern as well as several other rarely-measured parameters, such as photosynthetically active radiation and ratios of aboveground/belowground to total NPP.Thus, our method can be a more cost-effective tool for any practical efforts in managing, protecting, and/or restoring steppe grasslands, the very important but most beleaguered ecosystems on the planet [1].
The method was applied to the typical-steppe Balagaer River watershed located in northeast China.The results indicate that the evapotranspiration in this arid and semiarid watershed might be enhanced by grass canopy (αveg = 2.5 > 1) and compensated by soil water (αsm = 1.15 > 1).This is evidenced by previous studies (e.g., [9,28,55,56] showing that the annual actual evapotranspiration tends to be greater than the annual precipitation most of the years.Precipitation is probably not the sole water source for evapotranspiration.More soil water can be evaporated and transpired on the top of precipitation in a drier than a wetter year.A grass canopy can intercept portion of precipitation, providing more water available for evaporation, while grass roots can uptake soil water, making more water available for transpiration [41].Overall, the grazing removal rate in the area downstream of the West Ujimqin banner was higher, but had a smaller temporal variation, than that in the area upstream of the banner (Figure 5).The higher removal rate in the downstream area can be attributed to the larger population and thus more grass-grazing animals, whereas the lower removal rate in the upstream area is likely due to the lower population and thus fewer grass-grazing animals (Field reconnaissance by the authors in July 2014; [15]).Another possible reason is that the soils in the upstream area are less saline and more permeable than those in the downstream area,

Discussion
The method (Equations (1) through (11)) can be parameterized using easily-available data on topography, soil, land use and land cover, and climate.For an area of interest, the method predicts the natural NPP, whereas the MODIS images present the actual NPP.The removal rate can be estimated as the difference of the natural NPP less the actual NPP.In comparison to the conventional survey-based approach, the method is universally applicable and does not require animal-related data (e.g., species and heads).The method can be applied to areas with a small (or large) size and an unknown simple (or complex) grazing pattern.Leriche et al. [53] proposed a model for evaluating the short-term effect of grazing on grass NPP.This model requires inputs of grazing intensity and pattern as well as several other rarely-measured parameters, such as photosynthetically active radiation and ratios of aboveground/belowground to total NPP.Thus, our method can be a more cost-effective tool for any practical efforts in managing, protecting, and/or restoring steppe grasslands, the very important but most beleaguered ecosystems on the planet [1].
The method was applied to the typical-steppe Balagaer River watershed located in northeast China.The results indicate that the evapotranspiration in this arid and semiarid watershed might be enhanced by grass canopy (α veg = 2.5 > 1) and compensated by soil water (α sm = 1.15 > 1).This is evidenced by previous studies (e.g., [9,28,55,56] showing that the annual actual evapotranspiration tends to be greater than the annual precipitation most of the years.Precipitation is probably not the sole water source for evapotranspiration.More soil water can be evaporated and transpired on the top of precipitation in a drier than a wetter year.A grass canopy can intercept portion of precipitation, providing more water available for evaporation, while grass roots can uptake soil water, making more water available for transpiration [41].Overall, the grazing removal rate in the area downstream of the West Ujimqin banner was higher, but had a smaller temporal variation, than that in the area upstream of the banner (Figure 5).The higher removal rate in the downstream area can be attributed to the larger population and thus more grass-grazing animals, whereas the lower removal rate in the upstream area is likely due to the lower population and thus fewer grass-grazing animals (Field reconnaissance by the authors in July 2014; [15]).Another possible reason is that the soils in the upstream area are less saline and more permeable than those in the downstream area, leading to a downstream-to-upstream increase pattern of natural grass production [15].Averaged across the eleven analysis years, the watershed as a whole had a removal rate of about 47.6%, which is within the global HANPP (human appropriation of net primary production) ranges of 11% to 63% reported by [57].Here, the only different is that for our study watershed, grazing was the dominant grass removal mechanism, whereas for some other regions (e.g., Sothern Asia and Western Europe), both grazing and harvesting could be the grass removal mechanisms.For an area of interest, when pre-harvest MODIS images are available, our method may be used to better quantify grazing rate in the HANPP calculation, which is needed to quantify the aggregate impact of land use on biomass available each year in ecosystems.
In a previous study in the Balagaer River watershed, Wang et al. [42] conducted a sensitivity analysis of soil erosion to LAI.The take-home result from that analysis was reproduced and is shown in Figure 7 for the convenience of the discussion hereinafter.By assuming that a reference (i.e., an existing) LAI of 0.47 corresponds to the annual average watershed-scale removal rate of 47.6% and that the LAI for a given removal rate can be estimated by Equation ( 13), the fluvial erosion might not vary much but would be dominated by aeolian erosion (R aÑw > 1 in Figure 7b) during the analysis years.However, a removal rate of higher than 65% would drastically increase the magnitude and risk of soil erosion (steep increasing of R aw in Figure 7a), and thus can be treated as the threshold for grazing management in the Balager River watershed.

LAI "
1 ´premovalrateq 1 ´47.6% ˆ0.47 " 1 ´premovalrateq 0.524 ˆ0.47 " 0.8969 ˆr1 ´premovalrateqs (13) Water 2016, 8, 339 14 of 18 leading to a downstream-to-upstream increase pattern of natural grass production [15].Averaged across the eleven analysis years, the watershed as a whole had a removal rate of about 47.6%, which is within the global HANPP (human appropriation of net primary production) ranges of 11 to 63% reported by [57].Here, the only different is that for our study watershed, grazing was the dominant grass removal mechanism, whereas for some other regions (e.g., Sothern Asia and Western Europe), both grazing and harvesting could be the grass removal mechanisms.For an area of interest, when pre-harvest MODIS images are available, our method may be used to better quantify grazing rate in the HANPP calculation, which is needed to quantify the aggregate impact of land use on biomass available each year in ecosystems.
In a previous study in the Balagaer River watershed, Wang et al. [42] conducted a sensitivity analysis of soil erosion to LAI.The take-home result from that analysis was reproduced and is shown in Figure 7 for the convenience of the discussion hereinafter.By assuming that a reference (i.e., an existing) LAI of 0.47 corresponds to the annual average watershed-scale removal rate of 47.6% and that the LAI for a given removal rate can be estimated by Equation ( 13), the fluvial erosion might not vary much but would be dominated by aeolian erosion (Ra→w > 1 in Figure 7b) during the analysis years.However, a removal rate of higher than 65% would drastically increase the magnitude and risk of soil erosion (steep increasing of Raw in Figure 7a), and thus can be treated as the threshold for grazing management in the Balager River watershed.

Figure 1 .
Figure 1.Map showing the location, soils (i.e., World Food and Agriculture Organization or FAO soil classes), boundary, and subbasins of the Balagaer River watershed.

Figure 1 . 18 Figure 1 .
Figure 1.Map showing the location, soils (i.e., World Food and Agriculture Organization or FAO soil classes), boundary, and subbasins of the Balagaer River watershed.

Figure 3 .
Figure 3. Plot showing ratio of annual actual evapotranspiration (E) to annual precipitation (P), as well as E, versus P.

Figure 3 .
Figure 3. Plot showing ratio of annual actual evapotranspiration (E) to annual precipitation (P), as well as E, versus P.

Figure 4 .
Figure 4. Plots showing the: (a) MODIS-presented and model-predicted net primary production (NPP); (b) absolute NPP removal rate by grazing; and (c) percent removal rate, versus year.

2 )Figure 4 .
Figure 4. Plots showing the: (a) MODIS-presented and model-predicted net primary production (NPP); (b) absolute NPP removal rate by grazing; and (c) percent removal rate, versus year.

Figure 5 .
Figure 5. Plots showing the annual average and the minimum and maximum (error bars) of: (a) net primary production (NPP); and (b) NPP removal rate, for the study watershed as a whole and the 27 subbasins.Herein, the "banner" refers to the City of West Ujimqin, which is in the vicinity of the Baiyingwula flow station shown in Figure 1.

Figure 6 .
Figure 6.Plots showing net primary production (NPP) removal rate versus: (a) annual precipitation (P); (b) annual actual evapotranspiration (E); (c) E/P; and (d) atmospheric dryness (defined as the ratio of annual potential evapotranspiration E0 to P).

Figure 6 .
Figure 6.Plots showing net primary production (NPP) removal rate versus: (a) annual precipitation (P); (b) annual actual evapotranspiration (E); (c) E/P; and (d) atmospheric dryness (defined as the ratio of annual potential evapotranspiration E 0 to P).

Figure 7 .
Figure 7. Sensitivity of soil erosion to leaf area index (LAI) as influenced by soil moisture for: (a) total erosion Raw; and (b) aeolian in relative to fluvial erosion Ra→w.The definitions and formularizations of Raw and Ra→w can be found in [42].

Figure 7 .
Figure 7. Sensitivity of soil erosion to leaf area index (LAI) as influenced by soil moisture for: (a) total erosion R aw ; and (b) aeolian in relative to fluvial erosion R aÑw .The definitions and formularizations of R aw and R aÑw can be found in [42].

Table 1 .
The adopted values of the model parameters in this study † .

Table 3 .
The Pearson correlation coefficients between net primary production (NPP) removal rate and the major climate factors (P: annual precipitation; E: annual actual evapotranspiration; and E 0 : annual potential evapotranspiration).The bold coefficients are statistically significant at a significance level of α = 0.05.Notes: † See Figure1for the subbasins.UP: upstream of the City of West Ujimq; DS: downstream of the city.