The shortgrass steppe (SGS) occupies the southwestern part of the Great Plains. Half of the land is cultivated, but significant areas remain under natural vegetation. Despite previous studies of the SGS carbon cycle, not all aspects have been completely addressed, including gross productivity, ecosystem respiration, and ecophysiological parameters. Our analysis of 1998 — 2007 flux tower measurements at five Bowen ratio-energy balance (BREB) and three eddy covariance (EC) sites characterized seasonal and interannual variability of gross photosynthesis and ecosystem respiration. Identification of the nonrectangular hyperbolic equation for the diurnal CO2 exchange, with vapor pressure deficit (VPD) limitation and exponential temperature response, quantified quantum yield α, photosynthetic capacity Amax, and respiration rate rd with variation ranges (19 < α < 51 mmol mol-1, 0.48 < Amax < 2.1 mg CO2 m-2 s-1, 0.15 < rd < 0.49 mg CO2 m-2 s-1). Gross photosynthesis varied from 1 100 to 2 700 g CO2 m-2 yr-1, respiration from 900 to 3,000 g CO2 m-2 yr-1, and net ecosystem production from — 900 to 700 g CO2 m-2 yr-1, indicating that SGS may switch from a sink to a source depending on weather. Comparison of the 2004–2006 measurements at two BREB and two parallel EC flux towers located at comparable SGS sites showed moderately higher photosynthesis, lower respiration, and higher net production at the BREB than EC sites. However, the difference was not related only to methodologies, as the normalized difference vegetation index at the BREB sites was higher than at the EC sites. Overall magnitudes and seasonal patterns at the BREB and the EC sites during the 3-yr period were similar, with trajectories within the ± 1.5 standard deviation around the mean of the four sites and mostly reflecting the effects of meteorology.
Introduction
The shortgrass steppe (SGS) ecoregion occupies the southwestern part of the Great Plains of North America, covering approximately 0.34 106 km2 (Fig. 1). Half of the land is cultivated, but significant areas remain under natural vegetation dominated by the Grama-buffalo grass (Bouteloua-Buchloe) association (Lauenroth, 2008). As an important resource for agricultural production (cereals and animal products) and ecosystem services (Burke et al., 2008), the SGS ecoregion has been a focus of comprehensive systems — ecological studies, including the “Grassland Biome” project of the US IBP Program (Van Dyne, 1971). As a result, many aspects of the structure and functioning of the shortgrass steppe ecosystems have been thoroughly described (Lauenroth and Burke, 2008). In particular, biological productivity and element cycling of the shortgrass steppe have received special attention and eventually led to construction of dynamic ecosystem simulation models such as ELM and Century (Innis, 1978 and Parton et al., 1987, respectively). Several studies of CO2 exchange in ecosystems of the shortgrass ecoregion of North America and similar ecoregions of Europe and Asia were conducted using the chamber, flux tower, and remote sensing techniques (Brown and Trlica, 1977a; LeCain et al., 2000; Li et al., 2000; Gilmanov et al., 2004b, 2005, 2007, 2010; Fu et al., 2006, 2009; Belelli-Marchesini et al., 2007; Wu et al., 2008; Alfieri et al., 2009; Rey et al., 2012; Zhang et al., 2012; Rajan et al., 2013; Shao et al., 2013; Gao et al., 2014). Nevertheless, due to anticipated changes in climate and anthropogenic management, certain aspects of the SGS carbon cycle require additional scrutiny. Particularly, we do not have sufficient data on the ecosystem-scale estimates of fundamental characteristics of gross primary productivity (GPP), total ecosystem respiration (RE), and the resulting net ecosystem carbon budget, as the few available GPP and RE estimates of the North American SGS (Andrews et al., 1974; Brown and Trlica, 1977b; Detling, 1979; Risser et al., 1981; Klopatek and Risser, 1982) were based on extrapolating and modeling data from physiological studies at the leaf, plant, or chamber scales. Long-term measurements of ecosystem-scale CO2 exchange of SGS communities using the Bowen ratio-energy balance (BREB) and later the eddy covariance (EC) techniques began in the mid-1990s (Svejcar et al., 1997; Alfieri et al., 2009). There are different opinions concerning the evaluation and comparison of BREB and EC flux tower measurements of CO2 exchange of nonforest, particularly grassland ecosystems. Dugas et al. (1997) recognized BREB as an adequate tool for CO2-exchange measurements on grasslands, and the method was used by the US Department of Agriculture—Agricultural Research Service (USDA-ARS) Rangeland CO2 Flux project (Angell et al., 2001; Frank and Dugas, 2001; Sims and Bradford, 2001; Emmerich, 2003; Gilmanov et al., 2003b, 2006; Mielnick et al., 2005; Svejcar et al., 2008), as well as in other studies of grasslands and croplands (Dugas et al., 1991, 1999; Ham and Knapp, 1998; Asseng and Hsiao, 2000; Ansley et al., 2002; Gilmanov et al., 2003a; Baron et al., 2005; Scott et al., 2006; Irmak, 2010; Jamiyansharav et al., 2011; O'Dell et al., 2014). Phillips and Beeri (2008) have summarized long-term BREB measurements in the mixed grassland of North Dakota and established consistent relationships with remote sensing indices. Comparison of parallel BREB and EC system measurements demonstrated that although there are certain differences in energy and water fluxes, CO2 fluxes recorded by the two systems did not differ significantly (Dugas et al., 2001; Wolf et al., 2008). Hipps et al. (2002) observed reasonable agreement for water vapor fluxes measured by parallel running BREB and EC systems in a crested wheatgrass (Agropyron desertorum) field, wherein the EC CO2 fluxes were always larger than those measured with BREB. Nevertheless, results from both systems indicated that during the period of study the crested wheatgrass ecosystem was a net source of carbon. In contrast, Alfieri et al. (2009) found that BREB overestimates the magnitudes of carbon dioxide fluxes. Skinner and Wagner-Riddle summarized the problem: “Currently missing are studies comparing EC and BR flux estimates for the entire season or for complete annual cycles to determine how differences between systems affect long-term estimates of net C exchange” (Skinner and Wagner-Riddle, 2012, p. 377). Clearly, there is a need to compare seasonal patterns, annual totals, and ecophysiological parameters obtained from the two methods to evaluate opportunities to integrate the legacy BREB data accumulated from grassland and crops with the growing datasets from the EC networks.
Figure 1.
Major grassland ecoregions of the Great Plains (Omernik, 1987; Homer et al., 2015) and location of the study sites.

The objectives of this study, using all the available BREB and EC datasets of flux tower net CO2 exchange (F) measurements in ecosystems of the North American SGS ecoregion, are to 1) partition net CO2 fluxes into gross photosynthesis (Pg) and total ecosystem respiration (Re) components and estimate major ecophysiological parameters; 2) gap fill the data, describe seasonal patterns of CO2 exchange components and parameters, and estimate weekly and annual totals of gross primary production, total ecosystem respiration, and net ecosystem production (NEP); 3) compare CO2 exchange components and parameter estimates from BREB and EC flux towers, and 4) compare source/sink activity of ecosystems of the SGS ecoregion to mixed-grass and tallgrass ecoregions.
Materials and Methods
Study Sites
Five BREB and three EC tower sites considered in this paper (see Fig. 1) represent fundamental properties of grassland ecosystems of the SGS ecoregion and reflect features of the dominant management regimes: ungrazed, moderately grazed, and heavily grazed. Measurements at the BREB towers were conducted during 1998–2006 at several locations in northeastern Colorado and southeastern Wyoming. During 1998–1999, BREB tower measurements (location BR1, Table 1) were conducted at an ungrazed site on native shortgrass steppe at the Central Plains Experimental Range (CPER), administered by the USDA-ARS. In 2000 the BREB system was moved to a second ungrazed site at the CPER (location BR2). Both of these ungrazed sites were on an Olney fine-loamy soil (mesic Ustic Haplargids). The BREB station was moved a third time in 2004 to a heavy continuously grazed plot (BR3) located on a Remmit coarseloamy soil (mesic Ustic Haplocambids) at the CPER, and a second BREB tower was installed in a nearby moderate continuously grazed plot (location BR4) on a Zigweid fine-loamy soil (mesic Ustic Haplocambids). Characterized by a rather flat relief, the BREB sites have a spatially homogeneous vegetation of up to 0.3 m in height dominated by the C4 grasses blue grama (Bouteloua gracilis [H. B. K.] Lag. Ex Steud.) and buffalo grass (Buchloe dactyloides [Nutt.] Engelm.), accompanied by a mixture of C3 grasses (western wheat [Pascopyrum smithii {Rydb.} A.], needle-and-thread [Stipa comata Trin and Rupr.], and others); cacti; and shrubs. The grasses constitute > 70% of the total vegetation (Milchunas et al., 1989). According to both on-site measurement (LeCain et al., 2002) and remote sensing data (see corresponding section later), leaf area index at the site seldom exceeds 1.5 m2 m-2.
Table 1
Location of the flux towers in the SGS ecoregion of the North American Great Plains

The BREB site in southeastern Wyoming (BR5), located in the northwestern part of the SGS ecoregion, represents a true mixed-grass prairie dominated by the midheight cool-season grasses western wheatgrass and needle-and-thread grass. The site also contains a warm-season grass, blue grama, which is characteristic of shortgrass steppe. The soil of the site is an Ascalon sandy loam classified as a mixed, mesic, Aridic Argiustoll (LeCain et al., 2000). CO2 exchange measurements at BR5 were conducted during 1997–1998.
Ecological similarity of all sites is emphasized by the fact that in terms of the most detailed Level 4 Ecoregion taxonomy of the United States (Omernik and Griffith, 2008), sites BR1-BR5, EC1, and EC2 were classified as “High Plains. 25c. Moderate Relief Plains,” and site EC3 as “High Plains. 25i. Llano Estacado.”
The three sites using the EC methodology were within the native shortgrass steppe region but had been converted to seeded pastures or, at one time, row cropping. Measurements were conducted on ungrazed (location EC1) and moderately grazed (location EC2) USDA Conservation Reserve Program (CRP) pastures at the Curtis Ranch near Briggsdale, Colorado, from 2004–2007, and on a pasture seeded to warm-season C4 perennial bunch grass Caucasian bluestem (Bothriochloa bladhii [Retz] S.T. Blake) in the Texas High Plains (location EC3) from 2010–2011 (Rajan et al., 2013). Before 1987, the EC1 (and EC2) site had been in a wheat/fallow production rotation. In 1987, the site was placed in the CRP. It had no livestock grazing and a well-established cover of the cool season C3 grasses western wheat and needle-and-thread and the native warm season C4 grasses buffalo grass, sideoats grama (Bouteloua curtipendula [Michx.] Torr.), and little bluestem (Schizachyrium scoparium [Michx.] Nash). The soil of the site was classified as an Ascalon fine sandy loam (mixed, mesic, Aridic, Arjiustoll). In October 2003, a grazing treatment (EC2) was opened to spring and fall grazing at moderate intensity. Moderate grazing intensity was cattle grazing during two main periods in the spring and fall, until such time as approximately 250 kg/ha of forage remained in the pasture. Site EC3, located in the High Plains of Texas, was seeded to Caucasian bluestem in May 2007 and was grazed three times (May, July, August) in 2010, which was a high-productive year, but not grazed in 2011 due to extreme drought. The soil is a Pullman clay loam (fine, mixed, superactive, thermic Torrertic Paleustoll) with a flat relief (0%–1% slope) (Rajan et al., 2013).
Figure 2.
Mean annual temperature versus annual precipitation scatterplots (left) and precipitation histograms (right) for airport meteorological stations in (A) shortgrass steppe (95% confidence contour S), (B) mixed-prairie (contour M), and (C) tallgrass-prairie (contour T) ecoregions. Data were taken from the meteorological forecasting site Tutiempo Network S. L. Available at: http://en.tutiempo.net/climate/united-states.html.

The shortgrass steppe ecoregion has higher temperatures and lower precipitation (Fig. 2) than the mixed and tallgrass ecoregions, and most rains (∼70%) occur during May–September.
Meteorological conditions during the years of the study demonstrated a wide variety of weather patterns (Fig. 3). For example, at the ungrazed SGS Colorado sites (1998–2007), conditions varied from the wettest and coolest in 1999 (hydrologic year precipitation PCPNh = 545 mm, sum of temperatures above 5°C Tsum5 = 2 123-degree days) to the hottest and driest days in 2002 (PCPNh = 160 mm, Tsum5 = 2 344-degree days). At the Wyoming site, lower temperatures (MAT = 6.7°C) and higher precipitation (PCPNh = 437 mm) in 1997 resulted in higher net production (NEP = 436 g CO2 m-2 yr-1) than in 1998 (MAT = 7.8°C, PCPNh = 247 mm, NEP = 142 g CO2 m-2 yr-1). Weather conditions during the first (2010) year at the Lockney site were also marked by lower temperatures (MAT = 14.9°C) and high precipitation (PCPNh = 680 mm), followed by hot and dry conditions during the second (2011) year (MAT = 15.4°C, PCPNh = 187 mm). Such a wide range of meteorological factors provides a good opportunity to model climatic response of carbon dioxide exchange.
Instrumentation and Data Processing
Flux tower sites analyzed in this study were equipped with modern field instrumentation corresponding to the BREB or EC method described in Dugas (1993) and Campbell Scientific (1998) for BREB and Dugas et al. (2001), Meyers (2001), Alfieri et al. (2009), and Rajan et al. (2013) for EC. Standard correction procedures and outlier detection algorithms recommended for grassland ecosystems were applied to the raw datasets (Dugas et al., 1997, 2001; Falge et al., 2001 ; Wolf et al., 2008). In particular, during periods when the BREB algorithm was not valid for calculating turbulent diffusivity, it was estimated using atmospheric parameters as described by Dugas et al. (1999). Resulting “Level 2” (using Ameriflux terminology) files containing aggregated subhourly (20 min for BREB and 30 min for EC) values of the net CO2 fluxes (F) and the ancillary variables (incoming photosynthetically active radiation [Q], air temperature [Ta], soil temperature [5-cm depth, T s ], air relative humidity [RH], VPD, and others) served as inputs for the procedure of partitioning F into gross photosynthesis (Pg ) and ecosystem respiration (Re) components using the “light-soil temperature-VPD” response method (Gilmanov et al., 2013, 2014; Morgan et al., 2016). Following are the basic equations of the method:
Figure 3.
Seasonal and interannual dynamics of mean daily air temperature (T a ) and hydrologic year precipitation (PCPN h ) at the ungrazed shortgrass steppe site (A), Cheyenne site (B), and Lockney site (C).

where α is the initial slope (apparent quantum yield), Amax is the plateau (photosynthetic capacity) of the light-response, θ is the convexity parameter (Thornley and Johnson, 2000), r d is respiration rate when no temperature response was observed, r 0 and kT , are the coefficients of the exponential temperature response (r0 = Re[0]), and the normalized VPD-response function ϕ(VPD,VPDcr,σVPD) depends on two parameters: the critical value VPDcr, below which water deficit doesn't affect photosynthesis (ϕ= 1 for VPD ≤ VPDcr), and the VPD-response curvature parameter, σVPD (1 ≤ σVPD ≤ 30), with lower values describing a strong water-stress effect and higher values describing a weak effect (Gilmanov et al., 2013). The average daytime respiration rate rd was calculated as:
where t1 and t2 denote moments of sunrise and sunset, correspondingly. In addition, for days with identifiable Eqs. [1–4] parameters, exponential Eq. [3] was also fitted for the n-day window of the nighttime subhourly data centered at the day under consideration (n ≈ 9).
Figure 4.
CO2 uptake of the ungrazed shortgrass steppe ecosystem of the BR2 site for 2 June, 2003 (DOY 153): (A) light-response with the 95% confidence band (Eq. [6], α=0.00082 mg CO2 µmol -1; Amax =0.898 mg CO2 m-2 s-1; θ=0.409; r d =0.158 mg CO2 m-2 s-1); (B) light-soil temperature-VPD response (Eqs. [1–5], α=0.00077 mg CO2 µmol-1; Amax=1.07 mg CO2 m-2 s-1; θ=0.05; r0 =0.134 mg CO2 m-2 s-1, kT=0.009 °C-1, VPDcr=1.0 kPa; σVPD=4.3 kPa; gray dots—measurement data; black dots—model data; surface—response with mean daily VPD); (C) exponential regression of the nighttime respiration on soil temperature for DOY 149–157 with the 95% confidence band, r night = 0.074 exp(0.05 Ts ).

Parameters of Eqs. [1–5] characterizing diurnal dynamics of the CO2 exchange, such as apparent quantum yield α, photosynthetic capacity Amax, convexity of the light-response θ, and others, were numerically fitted to the datasets of individual measurement days of each measurement site-year. Interpolation and extrapolation of the seasonal patterns demonstrated by these parameters to days with missing measurements were used as major tools of gap filling (in addition to statistical interpolation of missing data for short subhourly intervals). Diurnal rates (mg CO2 m-2 s-1) of Pg , Re, and F were calculated by Eqs. [1–5] using the diurnal data for meteorological drivers (Q, Ta , Ts, RH, VPD). The 24-hr integration of these rates provided the year-round daily (g CO2 m-2 d-1) series of Pg(t), Re(t), and F(t) values (t = 1, 2, ..., 365) for corresponding years of study. Daily estimates of the ecosystem-scale ecophysiological parameters of photosynthesis and respiration in Eq. [1–5] were also obtained. For days when identification of Eq. [1–2] parameters was not possible (mostly outside the growing season), the net CO2 exchange for the day j was described as F(Ts) = -Re(Ts) with parameters of Eq. [3] estimated from the subhourly (F, Ts) data for the n-day window centered in day j (depending on data availability, n varied from 9 to 14).
Light-Use Efficiency
Among the diversity of light-use efficiency (LUE) coefficients (cf. Bonhomme, 2000), the two most frequently used are the physiological LUE coefficient calculated as a ratio of gross primary production P g to absorbed photosynthetically active radiation Qa , εphys = Pg/Qa, and the ecological LUE coefficient, εecol = Pg/Q, where Q is incident photosynthetically active radiation (Qa < Q). In contrast to εphys characterizing mostly physiological abilities of plants to assimilate atmospheric CO2, εecol also incorporates ecologically significant information about structure and architecture of the plant canopy, as εecol = εphys · fPARa(LAI), where fPARa(LAI) = Qa/Q is the fraction of PAR absorbed by the plant canopy. While physiological LUE coefficient εphys is considered a valuable characteristic providing a basis for the rapidly growing wave of publications on “production efficiency models” (PEMs) pioneered by Monteith (1972) and Sellers (1987), it was demonstrated that εphys poorly (Garbulsky et al., 2010) or even negatively (Runyon et al., 1994; Polley et al., 2011) correlates with LAI, evapotranspiration, and ecosystem productivity. In contrast, since the early period of plant production studies, εecol is known as a positive correlate of ecosystem productivity (Odum, 1959; Ničiporovič, 1968; Duvigneaud, 1974; Runyon et al., 1994). A particularly close relationship between ecological light-use efficiency and productivity was demonstrated for grasslands: Data by Sims and Singh (1978) for n = 36 site-yr of production measurements in grasslands of the western United States show highly significant correlation (r = 0.91, P < 10-6) between net primary production and ecological LUE. On these ecosystem comparisons, we are using the ecological LUE coefficient for brevity denoted below as ε = Pg/Q based on the incident PAR, Q.
Figure 5.
CO2 uptake of the ungrazed shortgrass steppe ecosystem of the BR2 site for 1 July, 2003 (DOY 182): (A) light-response with the 95% confidence band (Eq. [6], α = 0.00077 mg CO2 µmol-1; Amax = 0.193 mg CO2 m-2 s-1; θ = 1.0; rd = 0.115 mg CO2 m-2 s-1); (B) light-soil temperature-VPD response (Eqs. [1–5], α = 0.00082 mg CO2 µmol-1; Amax = 0.332 mg CO2 m-2 s-1; θ = 1.0; r0 = 0.254 mg CO2 m-2 s-1, kT = -0.041 °C-1, VPDcr = 1.0 kPa; σVPD = 3.16 kPa; gray dots—measurement data; black dots—model data; surface—response with mean daily VPD); (C) exponential regression of the nighttime respiration on soil temperature for DOY 178–186 with the 95% confidence band, rnight = 0.364 exp(-0.058 Ts ).

Remote Sensing Data
Remote sensing indices were recognized as useful tools for interpretation and scaling-up of on-site ecosystem-scale measurements (Wylie et al., 2004, 2007; Gilmanov et al., 2005; Heinsch et al., 2006). In this paper we used the 250-m data of the normalized difference vegetation index (NDVI) derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor for all sites and the 1-km resolution MODIS LAI estimates (DAAC/ORNL, 2015) as supplemental tools for comparing productivity of the BREE and EC sites in Colorado. More specifically, for the 2004–2006 period, we used the 7-day NDVI composites from the expedited MODIS (eMODIS) product (Jenkerson et al., 2010). The 8-d MODIS LAI values were recalibrated to match the 7-d eMODIS NDVI time step.
Results and Discussion
Light-Response Functions
Variability of the meteorological drivers affecting the SGS ecosystems during the measurement period (see Fig. 3) is translated into variability of their CO2 exchange. Two major environmental response patterns of the CO2 exchange may be distinguished as illustrated in Figures 4 and 5. The case of photosynthetic response dominated by incoming radiation Q is illustrated in Figure 4 by the BR2 site data for 2 June, 2003, which was a clear, sunny day (daily PAR total 51.9 mol m-2 d-1, Qmax = 2038 µmol m-2 s-1) with moderate temperatures (mean Ta = 14.2°C, Ta.max = 21.6°C) and low evaporative demand (VPDavg = 0.95 kPa, VPDmax = 1.68 kPa). As shown in Figure 4A, on this day, CO2 uptake of the SGS community is fairly well described by the univariate nonrectangular light-response function:
where rd is average daily ecosystem respiration rate and other parameters as described above in Eqs. [1–5]. This equation describes a significant part of the variance of the CO2 exchange (R2 = 0.95) achieving the mean squared error of SE = 0.063 mg CO2 m-2 s-1. There is little difference between the morning and afternoon branches of the light-response curve for day of the year (DOY) 153. Predominance of radiation as a major driver of CO2 uptake in this case is confirmed by the fact that switching from a univariate light-dependent function (Eq. [6], Fig. 4A) to the multivariate light-soil temperature-VPD response function (Eqs. [1–4], Fig. 4B) results in only a small additional improvement of the data fit (SE = 0.060 mg CO2 m-2 s-1, R2 = 0.96). Visually, the low significance of the VPD as a factor controlling the CO2 uptake in this case is reflected by the fact that the black dots in Fig. 4B, displaying model-generated F(Q, Ts,VPD) values, remain close to the response surface F(Q, Ts; VPDavg) (shown by a mesh) corresponding to the model [1]–[4] applied with a constant value VPDavg = 0.95 kPa equal to the average daily VPD for DOY 153. This surface, in its turn, is close to the gray dots indicating original data.
In contrast to day 153, day 182 (1 July, 2003) was a sunny, hot summer day (daily PAR total 53.68 mol m-2 d-1, Qmax = 2071.9 µmol m-2 s-1) with high temperatures (mean Ta = 20.1 °C, Ta.max = 35.6°C) and strong evaporative demand (VPDavg = 2.5 kPa, VPDmax = 5.2 kPa). The diurnal F(Q) plot for this day exhibited a strong hysteresis-like pattern with the afternoon branch significantly lower than the morning branch (see Fig. 5A). Temperature response of the nighttime ecosystem respiration rnight for the 9-d window centered on DOY 182 showed a decrease of respiration with negative exponential temperature response coefficient kT = -0.058 (°C)-1 (see Fig. 5C). This suggests that the decrease of the CO2 uptake in the afternoon likely cannot be explained by an increase of respiration with temperature, making vapor pressure deficit the most significant factor controlling F under drought stress conditions.
Figure 6.
Seasonal and interannual variation of the ecosystem-scale parameters at ungrazed shortgrass steppe (A), Cheyenne mixed prairie (B), and Lockney pasture (C): (a)—apparent quantum yield, α; (b)—light-use efficiency, ε; (c)—photosynthetic capacity, Amax ; (d)—daytime respiration rate, r d . Dots indicate weekly means; bars, errors of the means.

Table 2
Ecosystem-scale ecophysiological parameters of the SGS ecoregion ecosystems resulting from flux tower data processing by the light-soil temperature—vapor pressure deficit response method, Eqs. [1]–[6]

This conclusion is strongly supported by results of the nonlinear regression, which showed a highly significant (P < 0.0001) magnitude σVPD = 3.16 kPa of the curvature parameter in Eq. [2], at the same time indicating a negative value exponential respiration coefficient kT = -0.041 (°C)-1 of low significance (P value 0.076). VPD-response function [2] with parameters VPDcr = 1.0 kPa and σVPD = 3.16 kPa indicates a nearly fivefold decrease of CO2 uptake rate as the VPD on 1 July, 2003, increased from near 0 in the morning to 5.2 kPa in the afternoon.
Illustrating the VPD control of CO2 exchange on DOY 182, 2003, the black dots in Fig. 5B (full model [1]–[4] predictions) markedly deviate from the VPD-independent surface FVPDavg(Q, Ts) = F(Q, Ts, VPDavg) generated from Eq. [2], with the VPD fixed at its mean value for day 182, VPDavg = 2.5 kPa, and are close to the gray dots (measurement data).
Table 3
Maximum gross ecological light-use efficiency εmax (mmol CO2 mol incident quanta-1) of various types of semiarid grasslands

While plots similar to Figs. 4–5 were available for most of the days with valid {(Q, Ts , VPD, F)} datasets at a subhourly time step, a better way to demonstrate applicability of the [1]–[6] modeling scheme is to examine and evaluate the seasonal and interannual dynamics of the parameters resulting from identification of these equations (Fig. 6).
Seasonal and Interannual Dynamics of Ecophysiological Parameters
Time plots of the weekly aggregated ecophysiological parameters of the SGS ecosystems (Fig. 6) revealed both seasonal patterns of their change within the annual cycle and significant year-to-year variability.Parameters decrease by 50% to 67% in years with unfavorable conditions (yrs 2000, 2002, 2006, 2011, when either extreme drought occurred, PCPN < 180 mm, or low precipitation PCPN < 255 mm was combined MAT > 8.9°C) compared with years with higher precipitation and less heat stress (1998, 1999, 2001, 2003, 2010). Table 2 has a summary of parameter estimates of the SGS ecoregion resulting from flux tower data processing by the light-soil temperature-VPD response method. Depending on the meteorological and phenological conditions of the year, the seasonal curves of parameters may be either unimodal or bimodal (see Fig. 6). In SGS proper sites, maxima of up to 1.63 mg CO2 m-2 s-1 for daily photosynthetic capacity Amax and 45. 1 mmol mol-1 for the apparent daily quantum yield αmax were observed. In the bimodal years the midseason values may be as low as 0.2–0.5 mg CO2 m-2 s-1 for A max and 5–20 mmol mol-1 for αmax. These numbers compare well with the results by LeCain et al. (2003), who used steady-state leaf gas exchange measurements to evaluate native plants enclosed in open-top chambers. For the dominant C4 species Bouteloua gracilis, they estimated maximum Amax values of 1.14 for early season, 0.56 for midseason, and 1.05 mg CO2 m-2 s-1 for late season. For the subdominant C3 species Pascopyrum smithii, they obtained Amax values of 1.44, 0.35, and 1.08 mg CO2 m-2 s-1 for early season, midseason, and late season. According to Polley et al. (2010) maximum weekly Amax for the SGS site near Nunn, Colorado, is 0.81 mg CO2 m-2 s-1, close to our mean maximum of weekly Amax of 0.85 mg CO2 m-2 s-1 in Table 2. Tower-based Amax estimates are also available for two semiarid grasslands in Inner Mongolia, China, with annual precipitation range comparable with shortgrass steppe but lower annual temperature. For the Xilin Gol site (PCPN = 350 mm, MAT = -1°C) maximum Amax parameter was estimated at 0.55 mg CO2 m-2 s-1 (Fu et al., 2009), whereas for warmer Duolun site (PCPN = 386 mm, MAT = 2.1 °C) Amax = 1.32 mg CO2 m-2 s-1 (Zhang et al., 2012).
Figure 7.
Scatter diagrams of (A) apparent quantum yield versus photosynthetic capacity (A max , α) and (B) light-use efficiency versus daytime respiration rate (r d , ε) for ecosystems of the tallgrass, mixed-grass, and shortgrass ecoregions. The ellipses show the 95% statistical confidence zones for corresponding parameter sets. Data for the mixed-grass and tallgrass prairies are from Gilmanov et al. (2017a, b).

Figure 8.
Seasonal and interannual dynamics of gross photosynthesis, ecosystem respiration, net ecosystem production, and precipitation at (A) the Colorado ungrazed shortgrass steppe (1998–2007), (B) mixed-prairie, Cheyenne (1997–1998), and (C) seeded pasture, Lockney (2010–2011).

Ecosystem-scale quantum yield parameters in Table 2 are also in agreement with literature data for comparable communities. For an SGS site near Nunn, Colorado, Polley et al. (2010) found αmax = 55 mmol mol-1. Taking into account a certain overestimation of the slope of the light response by the rectangular hyperbolic model used by these authors (cf. Marshall and Biscoe 1980), this estimate matches well our maximum αmax = 51 mmol mol-1 (see Table 2). For dominant SGS species Bouteloua gracilis, (C4) and subdominant Pascopyrum smithii (C3) maximum quantum yield values for the early-season, midseason, and late-season periods were estimated as 47, 46, 66 mmol mol-1 and 66, 51, 53 mmol mol-1, respectively (LeCain et al., 2003). As the community-level value of α is usually lower than that of individual dominant species, these estimates compare well with αmax = 51 mmol. Using eddy covariance measurements at two dry steppe sites in Inner Mongolia with climatic conditions within the shortgrass steppe range (see Fig. 2A), John et al. (2013) estimated maximum initial slopes of the ecosystem-scale light response as 32.7 mmol mol-1 (Doulun site) and 41.6 mmol mol-1 (Xilinhot site), which is practically within the ± 1 standard deviation (SD) interval around the mean of 31.2 mmol mol-1 of daily αmax values in Table 2.
Table 4
Summary characteristics of the annual CO2 exchange of the shortgrass steppe SGS ecoregion ecosystems during 1998—2007

The range of the ecological light-use efficiency ε in Table 2, calculated as a ratio of total daily photosynthesis to total daily PAR, is approximately half that of α (Figs. 6A–6B), apparently due to the combination of light-saturation and water stress limitations of photosynthesis. Comparative data in Table 3 show that maximum ecological light-use efficiency values, ε max , of the Great Plains shortgrass steppe ecosystems fall within the range of estimates for semiarid grasslands of North America, Europe, and Asia, with shortgrass ε max being higher than in ecosystems from cooler and drier climates but lower than in communities with temperate climates (see Table 3).
The range of respiration rates, 0 ≤ rd ≤ 0.5 mg CO2 m-2 s-1, is approximately one-third the range of Amax (Fig. 6C–6D). Pulses of intense CO2 evolution observed in certain years (2000, 2002, 2003, ungrazed SGS, Fig. 6D) indicate occasional high metabolic activity of the SGS biota, often due to precipitation events (Austin et al., 2004; Parton et al., 2011; Fan et al., 2012).
As expected, ecophysiological parameters for the Cheyenne mixed prairie are comparable (A max, ε) or higher (α, rd ) than observed at SGS proper sites (see Table 2). The highest values of parameters were achieved at the seeded Caucasian bluestem (C4) pasture near Lockney, where in the rain-abundant year 2010, αmax = 51 mmol mol-1, ε max = 28.7 mmol mol-1, and Amax = 2.10 mg CO2 m-2 s-1.
To identify the place of SGS ecoregion ecosystems within the general ecophysiological parametric space of the Great Plains grasslands, we compared the (A max, α) scatter diagram and the (r d. ε) scatter diagram of the shortgrass ecoregion with the same plots for the mixed-grass and tallgrass ecosystems (Fig. 7A–7B). The ellipses in Figure 7 delineate the 95% confidence zones of each of the parameter sets. At both two-dimensional plots, the shortgrass ellipse is substantially different from the tallgrass ellipse, but it overlaps with the mixed-grass ellipse. Considering distribution of the parameters in the entire four-dimensional parametric space {(α ε, Amax r d )}, it turns out that the centroids µ S = (31.1, 16.2, 1.08, 0.31), µ M = (37.8, 21.4, 1.14, 0.29), µ T = - (48.1, 34.0, 2.14, 0.44) of the parametric clusters of the shortgrass, mixed-grass, and tallgrass ecoregions, correspondingly, are significantly different according to the Mahalanobis criterion (Rao, 1965) (P = - 0.0057 < 0.01 for short — mixed-grass comparison and much lower for comparison of short-tall and mixed-tall clusters).
Figure 9.
Gross primary productivity-ecosystem respiration (GPP-RE) scatter diagram for shortgrass steppe ecoregion sites (A) compared with the mixed-grass and tallgrass ecoregion sites (B). Contours S, M, and T denote 95% confidence zones for (GPP, RE) pairs for sites in the shortgrass, mixed-grass, and tallgrass ecoregions, respectively. Data for mixed and tallgrass sites are from Gilmanov et al. (2017a, b).

Parameters α ε, Amax, rd in Figure 7 using Eqs. [1]–[5] are ecosystem-scale (= canopy-scale as opposed to leaf-scale) parameters, meaning that, e.g., α is calculated per unit of incoming (not absorbed) PAR, and Amax measures photosynthetic capacity per m2 of ground surface (not per m2 of leaf area, as in physiological photosynthetic capacity Amax,L ). While there are complex interrelations among ecophysiological characteristics of plant canopies aimed at optimization of the photosynthetic uptake (Terashima and Hikosaka, 1995; Wright et al., 2005), leaf area, community diversity, and soil fertility are some of the major factors explaining lower α and Amax values in SGS compared with tallgrass prairies, with intermediate values in mixed-prairie communities (see Fig. 7). In particular, as demonstrated by Thornley and France (2007), the canopy Amax is proportional (directly proportional in monocultures) to the leaf area index LAI and Amax,L, resulting in lower photosynthetic capacity of SGS canopies with LAI ≤ 1.5 m2 m-2 compared with tallgrass canopies with LAI ≥ 3 m2 m-2. Similarly, lower LAI values in the SGS entail lower absorption of solar radiation, explaining lower canopy-level quantum efficiency α. Lower species diversity (including fewer N-fixers) and lower soil fertility in SGS compared with tallgrass prairies leads to lower N content, which also contributes to lower values of α and Amax (Evans, 1989; Peterson et al., 1999; Lee et al., 2003).
Seasonal and Interannual Dynamics of Photosynthesis, Respiration, and Net Ecosystem Production
The time series of daily photosynthesis Pg(t), respiration Re(t), and net CO2 flux F(t) obtained by application of Eqs. [1]–[5] to the raw tower flux data and their gap filling (Fig. 8A) confirms strong variability of the functioning of the SGS ecosystems emphasized earlier by production ecologists based on biomass and chamber CO2 exchange measurements and by range scientists based on forage and animal production studies. For example, LeCain et al. (2002) observed a twofold difference in maximum daily rate of CO2 uptake F in favorable and unfavorable years in the SGS at CPER in Colorado, with a sixfold difference between the midsummer minimum and the annual maximum during the year with bimodal pattern of F. Field estimates of the aboveground net primary production (ANPP) for the SGS at the CPER for the 1939–2014 period showed a broad range from 14 g m-2 yr-1 in 2002 (year of extreme drought) to 175 g m-2 yr-1 in 2009 (PCPNyr = 437 mm yr-1), with a mean of 93 g m-2 yr-1 and a coefficient of variation of 31% (data compiled from Lauenroth and Sala, 1992; SGS-LTER, 2013; Lauenroth, 2013; and Parton 2015, personal communication). The grazing season gain of yearling steers in the same ecosystem during 19 yr of study (1995–2003) was found to vary 2.4 times from 72 to 172 kg/head/season (Hart and Derner, 2008). Data for the Cheyenne (Fig. 8B) and Lockney (Fig. 8C) sites demonstrate the same pattern of strong dependence of grassland CO2 exchange on the precipitation regime and are in agreement with the observation by Derner and Hart (2007), who reported a fivefold difference in the peak standing crop in years with dry and wet springs in the mixed-grass prairie near Cheyenne, Wyoming, during 1991–2006.
Source-Sink Activity of Ecosystems of the SGS Ecoregion
Table 4 and Figure 9A show summary characteristics of the CO2 exchange in ecosystems of the SGS ecoregion. Annual photosynthetic uptake varied greatly from < 1 000 to nearly 4 000 g CO2 m-2 yr-1, driven mostly by moisture availability. At the same time, respiration losses varied from 900 to 3 700 g CO2 m-2 yr-1, resulting in performance variation from sinks with NEP up to 700 g CO2 m-2 yr-1, to sources with NEP as low as 900 g CO2 m-2 yr-1. As seen in Fig. 9A, the 95% confidence zone for the (GPP, RE) points of the SGS ecoregion is evenly bisected by the 1:1 diagonal. For comparison, a similar zone for the mixed-grass sites (Fig. 9B, contour M) has a larger portion of the confidence zone below the diagonal, indicating higher sink activity. The greatest proportion of data points falling below the diagonal occurs for the tallgrass sites (Fig. 9B, contour T), emphasizing a predominance of the preburn CO2-sink performance of the tallgrass prairies and the potential for significant accumulation of soil organic matter (Derner and Schuman, 2007).
Effects of Grazing
Statistical comparison between ungrazed and grazed site-years using the entire dataset (see Table 2) shows no significant difference (two-sided P values > 0.05) between means of gross primary production, ecosystem respiration, net ecosystem production, and eco physiological parameters Amax, rd.max αmax, and εmax (Fig. 10). The result remains unchanged (P values for differences of the means > 0.05) after exclusion of the outlier (site EC3 grazed in 2010 with high annual PCPN = 690 mm yr-1), except for annual respiration RE for which the difference between means of ungrazed and grazed site-years becomes significant (two-sided P value = 0.02 < 0.05). This result is in agreement with the observation by Milchunas et al. (2008) that 50-yr-long studies of grazing intensity treatments at SGS sites at the Central Plains Experimental Range (Colorado) showed average forage production rates of 75, 71, 68, and 57 g DM m-2 yr-1 for ungrazed, lightly grazed, moderately grazed, and heavily grazed treatments, respectively, range productivity being most sensitive to precipitation and soil fertility and, only last, grazing intensity. A comprehensive analysis of the effects of grazing and weather presented in Morgan et al. (2016) demonstrated that weather affects CO2 fluxes more than grazing practices in SGS ecosystems.
BREB-EC Comparison
The SGS ecoregion dataset provides an opportunity to compare BREB and EC systems data by either examining the concurrent flux measurements or comparing models and parameters on the basis of tower data from the two methods. Measurements from BREB and EC towers parallel and independently collected at similar (37-km distance) heavily and moderately grazed BR3 and BR4 sites and the ungrazed and moderately grazed EC1 and EC2 sites from 2004 to 2006, postprocessed by the same algorithm, [1]–[5], allow a comparison of the magnitudes and flux patterns at all four sites (Fig. 11). The curves of gross photosynthesis Pg for BREB sites are consistently higher than the curves for EC sites, while the respiration curve Re for the BREB sites is lower than for the EC sites. Corresponding 3-yr GPP cumulates for the BREB sites (4 106 and 4 588 g CO2 m-2) are higher than for the EC sites (3 676 and 3 982 g CO2 m-2), while the 3-yr RE cumulates for the BREB sites (3 792 and 3 789 g CO2 m-2) are lower than for the EC sites (4 970 and 4 920 g CO2 m-2) (Table 5). This might be considered as confirmation observations by Alfieri et al. (2009) that the BREB method generates higher net CO2 fluxes than the EC method. However, NDVI curves for the BREB sites lay mostly above the ones for the EC sites (particularly during the periods of maximum NDVI), emphasizing higher general greenness of the BREB compared with the EC sites (Fig. 11C). Comparing the seasonal NDVI integrals showed that 3-yr iNDVI totals for the BREB sites during 2004–2006 are significantly higher than for the EC sites for the same years (76.5 and 85.2 NDVI days vs. 61.3 and 45.2 NDVI d, P < 0.004 for all BREB-EC pairs). The LAI data (1-km resolution) demonstrate rather close dynamics at both tower types (Fig. 11D), with the 3-yr totals of leaf area duration (310 and 308 d) at the BREB sites somewhat higher than at the EC sites (309 and 300 d).
Figure 11.
Seasonal and interannual dynamics of the weekly photosynthetic uptake (A), ecosystem respiration (B), eMODIS normalized difference vegetation index (C), and leaf area index (D) at the high continuously grazed (BR3), moderate continuously grazed (BR4), ungrazed (EC1), and moderately grazed (EC2) shortgrass steppe sites during 2004–2006. “Mean” indicates mean of the four sites (solid line); “Low,” mean is 1.5 standard deviation; “High,” mean is + 1.5 standard deviation.

Table 5
3-yr integrals of the normalized difference vegetation index (NDVI) truncated at the background, leaf area duration, gross primary production, ecosystem respiration, and net ecosystem production at the Bowen ratio-energy balance and eddy covariance sites, 2004–2006.

Comparison of ecophysiological parameters estimated from BREB and EC towers (Fig. 12) showed no significant differences between centers of the scatter ellipsoids of the BREB-derived and EC-derived parametric clusters.
In the entire four-dimensional space of maximum parameter values {(α ε Amax, r d )}, Mahalanobis criterion (Rao, 1965) showed no statistical grounds to reject the hypothesis that the BREB-based centroid µ BREB = (31.9, 16.3, 1.11, 0.33,) is not significantly different from the EC-based centroid µ EC = (29.8, 15.9, 1.03, 0.26) (P = 0.243 » 0.05). In other words, there is no significant difference between parameter means obtained from BREB or EC towers.
Thus, our study confirms the conclusion by Wolf et al. (2008) that the results of flux measurements by the BREB and EC towers are not essentially different, though “subtle differences” do occur under certain conditions. Figures 11A and 11B show that during the growing season, both of the Pg(t) and Re(t) curves remain within the {±1.5 SD} band around the mean of the combined BREB and EC data. The EC respiration estimates outside the growing period are definitely higher than the BREB-based estimates (see Fig. 11B), contributing to the overall higher RE totals for EC compared with BREB (see Table 5). Both methods have technical difficulties during the cold season (Dugas et al., 1999; Gilmanov et al., 2004a; Burba et al., 2008). The higher gross and net CO2 uptakes recorded at the BREB sites compared with the EC sites might be related to not only the higher respiration losses at the EC sites but also possibly higher general productivity of the BREB sites (BR3 and BR4), as indicated by their higher NDVI and LAI indices (Fig. 11C and D) and higher integrated NDVI and leaf area duration (see Table 5).
Summarizing the rates and patterns of the CO2 exchange in rangelands of the SGS ecoregion, this study demonstrated considerable variability of daily maxima for photosynthesis 13 < Pg < 41 g CO2 m-2 d-1, for respiration 13 < Re < 35 g CO2 m-2 d-1, and for net daily CO2 flux 9 < F < 24 g CO2 m-2 d-1. The annual ranges were 1 100 < GPP < 2 700 g CO2 m-2 yr-1, 900 < RE < 3 000 g CO2 m-2 yr-1, and -900 < NEP < 700 g CO2 m-2 yr-1. Both the daily and the annual fluxes in shortgrass steppe were lower than in mixed-grass and tallgrass ecosystems. Depending on meteorological conditions each year (chiefly precipitation amount and distribution), either unimodal or bimodal patterns of the seasonal CO2 exchange were observed and the ecosystem response varied from a net sink of + 700 g CO2 m-2 yr-1 in 1999 with 487 mm precipitation to a net source of -900 g CO2 m-2 yr-1 in 2010 with 173-mm precipitation. Parameterization of CO2 exchange using the nonrectangular hyperbolic equation with VPD limitation and exponential dependence of respiration on soil temperature provided quantitative estimation of apparent quantum yield α, ecological light-use efficiency ε, photosynthetic capacity Amax, and ecosystem respiration rate r d , including determination of the ranges of variation (19 < α < 51 mmol mol-1, 7 < ε < 29 mmol mol-1, 0.48 < Amax < 2.1 mg CO2 m-2 s-1, 0.15 < r d < 0.49 mg CO2 m-2 s-1) and patterns of their seasonal and interannual dynamics. Numerical values of parameters derived from flux tower measurements agree with estimates obtained by other authors using the open-top chamber technique. While there were certain differences in CO2 fluxes between BREB and EC towers, they occurred occasionally and were not related only to different methodologies because the NDVI and LAI values at the BREB sites were generally higher than at the EC sites. These differences were not significant in terms of the overall seasonal patterns (both systems responded to meteorological drivers with trajectories of Pg and Re remaining within a ± 1.5 SD band around the mean) or in the magnitudes of annual GPP, RE, and NEP totals and parameters. This suggests a need to consider including the legacy of the BREB CO2 exchange data, especially extensive in rangeland ecosystems, in ongoing efforts of comparative analysis, synthesis, and upscaling of rangeland CO2 exchange and productivity.
Figure 12.
Comparison of parameters estimated from Bowen ratio-energy balance and eddy covariance (EC) towers: (A) maximum apparent quantum yield α versus maximum photosynthetic capacity Amax ; (B) maximum light-use efficiency ε versus maximum daytime respiration rate rd . Blade quadrat indicates mean of the BREB estimates; black circle, mean of the EC estimates; thick ellipses, 95% zone for BREB parameters; thin ellipses, 95% zone for EC parameters.

Implications
The light-soil temperature-VPD-based method provides consistent partitioning of the raw flux tower CO2 exchange measurements in shortgrass rangelands into gross primary production (GPP) and total ecosystem respiration (RE) components while avoiding overestimation of respiration inherent to the widely used partitioning method based on nighttime temperature response functions. No significant differences were observed between GPP and RE estimates and CO2 exchange parameters at SGS sites with BREB and EC tower types. Long-term flux tower measurements demonstrate that SGS ecosystems switch from a sink to a source for atmospheric CO2 depending on weather. Quantification of the rates and parameters of the SGS ecosystems presented in this paper strengthens the empirical basis for sBS0patiotemporal modeling and upscaling as tools for forage production and carbon management of the SGS rangelands.
Acknowledgments
The authors thank Dr. Justin D. Derner for consultation and managing the heavy and moderately grazed pastures at the CPER, Dr. Daniel R. LeCain for assistance with BREB and plant biomass measurements, and Dr. Willam J. Parton for consultations in quantitative ecology of SGS ecosystems.
References
Notes
[1] ☆The study was supported in part by the USDA-ARS global change program for carrying out research; SGS-LTER for funding some of the data processing, analysis, and write-up; NIGEC-DOE project Carbon, water and land-use in Conservation Reserve Program lands of the shortgrass prairie for measuring and data processing of the Colorado EC flux tower sites; USGS Land Change Science Program, USGS Contract G15PC00012 for support of remote sensing data processing.