Before European settlement, the Great Plains of the United States contained vast herds of bison. These bison altered the landscape through their grazing. Measurement data of the disturbance that such grazing could produce, when scaled for the large population of bison, were used with a coupled atmospheric–ecosystem model to evaluate the likely effect that this grazing had on the growing season weather in the Great Plains. A dynamically coupled meteorological and plant growth model was used to investigate the regional atmospheric conditions over a single growing season. A 50-km horizontal mesh was implemented, covering the central plains of the United States. The modeling system was then integrated, with a time step of 90 s, for a period covering 1 April 1989 through 31 August 1989 using boundary conditions obtained from an objective analysis of gridded archive data. This integration was performed with and without grazing to assess the effects on regional atmospheric and biological processes. The grazing algorithm was employed to represent presettlement North American bison and was switched on and off for different simulations. The results indicated a cooling response in daily maximum temperatures to removal of grazing. The opposite trends were found for the minimum daily temperature. It was also found that grazing produced significant perturbations in the hydrological cycle.
Prior to the mid-1800s, tens of millions of bison lived on the Great Plains of the United States (Dang 1997). These bison grazed vast areas of this region. Given that it is well documented that grazing alters landscape (e.g., NRC 1992), an important research question is whether the change in above-ground biomass associated with grazing could have had a significant effect on the weather during the growing season. Measurements of actual grass consumption by individual bison using current-era data can be used to estimate the alteration of the Great Plains grasslands by the tens of millions of bison.
Over the last decade, progress in vegetation modeling suggests that simplified assumptions of static vegetation in weather and climate models are questionable (Pielke 1998; Neilson and Drapek 1998). The scientific community is now turning more attention to the role of land cover on regional and global climate (Chase et al. 1996, 1999; Eastman et al. 2001; Lu et al. 2001; Stohlgren et al. 1998). Chase et al. (1996), for example, found that regional changes in vegetation leaf area index (LAI, m2 m−2) affect global circulation patterns.
There is also ample evidence in recent modeling efforts that suggests a primary role of the biosphere when simulating paleoclimatology (Claussen and Gaylor 1997; Ganopolski et al. 1998; deNoblet et al. 1996; Crowley and Baum 1997; Texier et al. 1997; Brostrom et al. 1998). Most of these models simulate equilibrium conditions in which a year or more of atmospheric conditions are fed to an ecological model which is then allowed to come to equilibrium. They are termed asynchronous modeling systems. These models do not yet incorporate short-temporal-scale feedbacks and are often too coarse to resolve regional-scale spatial effects.
The current terrestrial biosphere is largely influenced by human interaction. This includes agricultural practices and grazing. In addition, the growth of agricultural crops occurs on relatively short timescales, often over a single growing season. To capture the short-term and regional response requires a linkage of physically consistent, fine-resolution plant and atmospheric models, as well as considerable computational resources. This study will address grazing effects at these finer spatial and temporal scales.
The effects of overgrazing have been seen in many ecosystems worldwide. Areas such as the Sahel, Mongolian short grass steppe, and the North American Great Plains have experienced intense grazing pressures in the past, and grazing continues to be significant (NRC 1992). Because grazing pressure is analogous to a change in LAI such as was investigated by Chase et al. (1996), it could be expected that grazing could alter large-scale and regional atmospheric circulation patterns, which would then affect ecosystems and plant communities.
The common practice of using prescribed surface boundary conditions or equilibrium models is not adequate for capturing biological feedback on climate. In this paper we present a highly mechanistic, physically based modeling system that is capable of simulating the interactions between atmosphere and biosphere. We have linked the Regional Atmospheric Modeling System (RAMS; Pielke et al. 1992) with the General Energy and Mass Transfer Model (GEMTM; Chen and Coughenour 1994; Coughenour and Chen 1997), producing a model we now refer to as “GEMRAMS.” In this study we use GEMRAMS to examine a seasonal, regional response to grazing. An earlier paper (Eastman et al. 2001) used this modeling system to investigate the effects of the conversion of natural to current landscape and increased levels of carbon dioxide (CO2) on a growing season in the same region as this study.
2. Modeling system
The meteorological model used in this investigation is based on RAMS. For a complete description of RAMS, the reader is referred to Pielke et al. (1992). RAMS, the spatially explicit mesoscale model, was further modified (Liston and Pielke 2000) for longer integrations and is named “ClimRAMS.”
The ClimRAMS model is a three-dimensional primitive equation model. The model is typically integrated over monthly to yearly timescales with horizontal grid spacings on the order of 50 km covering continental scales. In contrast, a typical general circulation model uses grid spacing on the order of several hundreds of kilometers. ClimRAMS is also capable of integrations ranging from a meter or less to several thousand kilometers. The vertical grid spacing can also be adjusted to address the specific problem being examined. In this case, the vertical mesh was 125 m at the surface and was stretched to 1500 m above the tropopause. A 90-s time step is used for the entire integration. At every time step and grid point, the model simulates three-dimensional momentum, temperature, and moisture fields. Also updated at each time step are precipitation, soil moisture and temperature, turbulence, and longwave and shortwave radiation, including the effects of clouds.
The boundary layer parameterization is based on similarity theory. The plant model that is coupled to the meteorological model also uses similarity theory, which maintains consistency between the models. Turbulence is parameterized through an anisotropic deformation-type K scheme.
The precipitation is parameterized by the Kuo cumulus convection scheme, in which a minimum vertical motion at the diagnosed cloud base is used to trigger convective precipitation, and by a large-scale precipitation algorithm using a dump-bucket precipitation scheme. This algorithm diagnoses supersaturation within the cloud and removes excess water based on precipitation efficiency. The water removed by these parameterizations is then received at the surface, where it is used in the soil moisture and energy budget. If the soil is already at saturation, the water is considered to be runoff.
The land surface model in ClimRAMS is broken into four compartments. Separate surface energy budget calculations are performed for the water, bare and shaded soil, and the dominant vegetation contained with a grid cell. In this model, vegetation is modeled as a biome type. Figure 1 indicates the biome types used in this model. These separate calculations are then combined through linear weighting. The vegetation calculations are based on a single-layer canopy model (Avissar and Mahrer 1982). The stomatal conductance is a function of relative humidity, incoming solar, surface temperature, soil water, and ambient CO2. The calculation of stomatal conductance is done every time step with GEMTM, the plant model coupled to ClimRAMS.
GEMTM explicitly calculates C3 (Fahquhar et al. 1980) and C4 (Collatz et al. 1992; Chen and Coughenour 1994) photosynthesis rates during each time step. The photosynthesis rates are computed for diffuse and direct photosynthetically active radiation (PAR). PAR is computed using a multilevel canopy radiation model (Goudriaan 1977). Photosynthesis rates are also dependent on plant water potential (PWP), vegetation temperature, atmospheric CO2, and the mixing ratio of water vapor. Instead of using a function for each of these variables to calculate photosynthesis and subsequent stomatal conductance, a highly mechanistic algorithm is used to simulate their effects. For a complete description of these dependencies the reader is referred to Chen and Coughenour (1994). Also during the time step, allocations of photosynthate, respiration, death rate, and growth rate are updated based on air temperature, soil temperature, and soil moisture relations. For allocation, the water status of the roots is used to determine which fraction of a biome-specific potential allocation will be partitioned to the roots. The remaining photosynthate is then partitioned to the leaves, stems, and seed, based on the phenology of the plant. The allocation coefficients used are also applied to the growth rates. These algorithms have been successfully validated for several grasses (Chen et al. 1994). The maximum death and growth rates are modified through temperature and water status by interpolation. The model uses four input temperature and water status values, each with an associated percentage of maximum death rates. These are then multiplied together, along with the maximum death rate, to calculate the actual death rate.
Tightly coupled to the hydrological cycle is a spatially explicit root model for effluence and uptake (Chen and Leith 1992), including branching and lengthening algorithms. The root model is embedded with the ClimRAMS soil model (Tremback and Kessler 1985). The GEMTM root model simulates root-level processes in greater detail than the default root model in RAMS, which extracts moisture from the moistest layer in the soil. The GEMTM root model extracts/efflues soil moisture at all levels for which root biomass is present (effluence is commonly referred to as hydraulic lift, which has been observed in numerous plant species). The root resistance to moisture transport depends on the radial resistance to water flow across the root to the xylem, a resistance in the xylem, and a soil resistance that is dependent on the grid cell's soil type. The first two resistances are parameterized using root length density (RLD) as a dependent variable. RLD is updated daily, based on the allocation of photosynthate to the roots used in conjunction with a branching and lengthening algorithm. In this way the plants are allowed to forage for water. The soil resistance is also dependent on RLD and on the time-dependent soil hydraulic conductivity. The model makes the assumption that transpiration and soil water uptake are in balance, which allows a solution to the PWP and integration with photosynthesis calculations.
Last, net photosynthate is accumulated every time step and is used at the end of a 24-h period in determining updated LAI, fractional coverage, and roughness length, all of which are used in the next 24 h of integration. Thus LAI, fractional coverage, and roughness length are used to modify the surface heat and moisture flux.
3. Experimental design and initialization
The coupled modeling system, hereinafter referred to as GEMRAMS, was used for two simulations: with and without grazing. For all integrations, the horizontal extent of the modeling domain and the vegetation distribution are shown in Fig. 1. The grid consists of 30 points in the east–west direction, 36 north–south points, and 20 vertical levels. The horizontal grid increment was 50 km, and the vertical mesh was stretched from 200 m near the surface to 1 km near the model's top, which was at 21 km. All meteorological fields were initialized using objective analysis software to process the National Centers for Environmental Prediction–National Center for Atmospheric Research weather reanalysis data (Kalney et al. 1996; Pielke et al. 1998) for 1 April 1989. Subsequent boundary conditions were derived using the same software and were applied at the outer three grid points through a Newtonian relaxation algorithm at 12-h intervals from 1 April 1989 through 31 August 1989. These fields provide identical boundary conditions for all integrations, leaving only grazing as the forcing mechanism to produce the differences between the simulations. In this way, the investigator can determine the model's sensitivity to a given forcing without concern for meteorological differences advecting into the domain from the outside.
The model was validated for the temporal period simulated in this study by comparing its output with surface observations and LAI derived from the 8-km Advanced Very High Resolution Radiometer (AVHRR) data archived at the Earth Observing System data center. The validation simulation was performed using current land use patterns as the vegetation cover. A correlation coefficient of 0.89 was computed for the maximum daily temperature, with a bias of −0.98°C, as reported in Eastman et al. (2001). The results were slightly better for the minimum daily temperature, for which the correlation coefficient was 0.91 and the bias −0.16°C. The precipitation correlation was found to be only 0.35, with a bias of −0.60 mm. Although the model did a fair job of simulating large-scale synoptic precipitation, it underpredicted the total precipitation during summertime convective situations. This is believed to be an artifact of the coarse grid resolution (grid spacing of 50 km), weaknesses in the simple microphysics parameterizations, and difficulty in scaling up point observations of precipitation to the model grid. Also, the modeled LAI was compared with the AVHRR-derived LAI. It was found that the model correlation coefficient was 0.90, with a bias of −0.19 m2 m−2.
The AVHRR-derived LAI was also used to initialize the leaf, shoot, and root biomass fields. Separate algorithms were used for woody LAI (Nemani and Running 1989) and herbaceous LAI (Asrar et al. 1984), which were then converted to leaf biomass using a biome-unique specific leaf area index (m2 kg−1). The shoot and root biomass fields were then initialized using leaf-to-shoot and shoot-to-root ratios specified for each biome type. A maximum rooting depth and depth of maximum root biomass were then used to distribute vertically the root biomass through linear weighting.
The grazing algorithm is designed to mimic the level of grazing that would be imposed by the herds of bison that roamed the North American grasslands in presettlement times (Dang 1997). It was assumed that 40 000 000 buffalo (National Bison Association, 2000) carried an average weight of 450 kg per head and that their daily intake was 2% of their body weight (Hawley et al. 1981), or 9 kg. It is well known that bison moved from north to south during autumn then returned in the spring as grass began growing again, generally along circular routes covering around 300–600 km. In an attempt to mimic this behavior, grazing was distributed in proportion to LAI by weighting the LAI values (a minimum of 0.1 was required) at grassland grid cells every 10 days. Thus, in the early spring portion of the integration, when grass was still senescent in the northern parts of the domain, the bison were primarily in the southern portion of the domain. As the northern grasses continue to grow there is a corresponding movement to the north. Although the algorithm is relatively simple, it seems to simulate the general seasonal migration of bison.
4. Results and discussion
a. Domain averages
We first examined domain-average values for 150 days over all grid points except the nudged lateral boundaries. The averages were also smoothed temporally with a 5-day sliding window. This procedure was repeated for several meteorological and biological variables. The first, maximum daily temperature (°C) indicates a maximum cooling of nearly −0.7°C when grazing is not included, with a minimum of 0.0°C at the start of the integration. The cooling effect increases steadily during the integration. This would suggest that as the bison remove biomass throughout the integration the transpiration is reduced. The reduction in transpiration would then result in a decrease in the latent heat surface turbulent fluxes and subsequent increases in sensible heat, which leads to higher daytime temperatures.
Examination of the domain averages for minimum daily temperature (Fig. 2b) indicates small warming trends, although the data show some oscillation around no contribution. At night the greater coverage of the vegetation when it is not grazed leads to higher atmospheric water vapor content, which elevates the dewpoint temperature, thus leading to warmer minimum temperatures.
Figures 2c and 2d display the effects on total accumulated precipitation and total accumulated convective precipitation, respectively. The total accumulated precipitation includes convective and large-scale precipitation. During the first month of integration we can see that the effects are minimal for both precipitation fields. As the growing season proceeds and summertime temperature increases, the effects are noticeable. Similarly, Eastman et al. (2001) found that a decrease in daily maximum temperature leads to a decrease in convective precipitation. The heating that lifts a parcel of air to the level of free convection is reduced at these lower temperatures when grazing is excluded.
Domain-average differences for biological data are shown in Figs. 3a–d. In Figs. 3a and 3b the data display a steady increase in leaf biomass when no grazing is present. The plot for root biomass also shows a similar trend in response to the grass trying to maintain a fixed shoot-to-root ratio. As the bison remove above-ground biomass, the plant compensates by allocating more photosynthate to replace the removed above-ground biomass. When grazing is absent, the grass does not have to compensate for the removal. Maintenance respiration, shown in Fig. 3d, indicates a higher value when grazing is not present, corresponding to the maintenance respiration's dependence on biomass. There is also a temperature effect in this calculation that would tend to increase the maintenance under warmer temperatures. This is the likely contributor to the sign reversal in the plot that occurs in early August.
b. Spatial analysis
Here, selected variables are averaged over the integration period and displayed as spatially explicit fields. Figures 4a and 4b, 5a and 5b, 6a and 6b, 7a and 7b, and 8a and 8b show the control simulation in the top panel and the effect of removing the bison in the bottom panel. The maximum daily temperature (°C) for the control simulation, displayed in Fig. 4a, indicates the topographical effects expected for this domain, with cooler temperatures at higher elevations. There is also a hint of a north–south gradient over the central plains. In Fig. 4b, there is a clear impact due to the lack of bison grazing over the plains. Maximum daily temperatures are 0.4°–0.8°C cooler over the plains, with a few small patches even cooler. There are also small changes of ±0.2°C in the mountainous regions to the west, although these perturbations are small. It was found in Eastman et al. (2001) that there was a biological teleconnection between those areas where land use had changed from a state of natural vegetation to current vegetation and those where there had been no change. Despite no land use change, these points had a different local climate. Here we can see points where there is no grazing, indicating a change in daily maximum temperature, albeit small. There is a gradient along the Front Range of the Rocky Mountains where the cooling is between 0.2° and 0.4°C when grazing is removed. This suggests that the teleconnection is a communication of atmospheric temperature and moisture fields through wind advection.
In the next series, Figs. 5a and 5b, the control simulation and grazing effects are shown for minimum daily temperature. The control simulation, Fig. 5a, exhibits spatial patterns similar to those seen for the maximum daily temperature, with an enhancement in the north–south temperature gradient over the plains. In the case of the no-grazing effect, Fig. 5b, the east–west boundary is not as apparent as in the maximum temperature plots. Overall, there is a mixture of cooling and warming throughout the domain; however, the magnitude of the warming is slightly greater. It was also found in the domain-averaged plots that the signal was somewhat mixed and oscillated between warming and cooling.
The analysis of domain averages suggests that the accumulated convective precipitation exhibited a large effect due to the lack of grazing. Figures 6a and 6b show the spatial distribution of precipitation in the control simulation and the effect on precipitation due to no grazing. The control distribution indicates an east–west gradient in convective precipitation. For the grazing effect in Fig. 6b, there is a reduction over the plains of between 5 and 20 mm. As suggested by the domain averages, it is likely that this reduction is due to the cooler daytime temperatures. There is also a general increase in the mountainous regions of between 0 and 5 mm. In these areas the dominant forcing mechanism for precipitation is topography.
In Fig. 7a (the leaf area index for the control simulation), the contrast between the plains and mountainous regions is clearly indicated by low and high LAI values, respectively. The no-grazing effect displayed in Fig. 7b shows a general increase in LAI of between 0 and 0.2 over the plains, and the LAI has changed by between ±0.2 in the mountainous regions. The increase over the plains is more significant, because the LAI for grasses is on the order of 1, while LAI in many of the mountainous biomes is on the order of 10. These slight increases and reductions in the western portion of the domain are close to the noise level of the model.
The domain and spatial analyses discussed in the previous section outline a complicated meteorological and biological response of biomes to grazing. The most significant result was the large cooling response in daily maximum temperatures due to the removal of grazing. It was also suggested that the cooler daytime temperatures could lead to reduced precipitation.
The minimum daily temperature fields exhibited a mixed trend. Overall, there was a warming in nighttime temperatures when grazing is removed, although there is variability in the domain averages. The increased water vapor content is a possible mechanism for the slight increase in minimum temperatures, an effect that would be minimized when there is strong synoptic-scale moisture transport.
Biological response was investigated for several related variables. Grazing impacts were important for above- and below-ground biomass. The root response is most intriguing, and its effect on the hydrological cycle needs further investigation. When grazing is present, the ground biomass is removed, which would lessen the demand on water resources. At the same time, the plant will attempt to compensate by allocating less photosynthate to roots, which reduces the root biomass and the ability to remove soil water. The effect of this over several years needs examination.
Overall, this study answers the question, “Does grazing affect regional climate?” with a “yes.” We conclude that the impacts of grazing are a significant component of the biosphere and should be accounted for in areas where grazing pressures are (or had been) large. However, the results of this study need to be taken in context. Whether these regional-scale impacts persist over decades or more is an important topic for further research that can be answered as computational efficiency increases.
Acknowledgments are extended to the institutions that provided the necessary support to complete this work. These include NASA Grant NAG8-1511, NSF Grant ATM-9306754, EPA Grant R824993-01-0, and NSF Grant DEB-9524129. Dallas Staley and Allison Edwards completed an excellent editorial improvement to the paper.
Current affiliation: Department of Atmospheric Science, Colorado State University, Fort Collins, Colorado.
Corresponding author address: Joseph L. Eastman, Dept. of Atmospheric Science, Foothills Campus, Colorado State University, Fort Collins, CO 80523. Email: email@example.com