Productivity and CO2 Exchange of Great Plains Ecoregions. I. Shortgrass Steppe: Flux Tower Estimates☆

ABSTRACT 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 10 6 km 2 (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 , 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 . 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 andParton et al., 1987, respectively). Several studies of CO 2 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., 2004bGilmanov et al., , 2005Gilmanov et al., , 2007Gilmanov et al., , 2010Fu et al., 2006Fu et al., , 2009Belelli-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 CO 2 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 CO 2 exchange  (Omernik, 1987;Homer et al., 2015) and location of the study sites. of nonforest, particularly grassland ecosystems. Dugas et al. (1997) recognized BREB as an adequate tool for CO 2 -exchange measurements on grasslands, and the method was used by the US Department of Agriculture−Agricultural Research Service (USDA-ARS) Rangeland CO 2 Flux project (Angell et al., 2001;Frank and Dugas, 2001;Sims and Bradford, 2001;Emmerich, 2003;Gilmanov et al., 2003bGilmanov et al., , 2006Mielnick et al., 2005;Svejcar et al., 2008), as well as in other studies of grasslands and croplands (Dugas et al., 1991(Dugas et al., , 1999Ham 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, CO 2 fluxes recorded by the two systems did not differ significantly 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 CO 2 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.
The objectives of this study, using all the available BREB and EC datasets of flux tower net CO 2 exchange (F) measurements in ecosystems of the North American SGS ecoregion, are to 1) partition net CO 2 fluxes into gross photosynthesis (P g ) and total ecosystem respiration (R e ) components and estimate major ecophysiological parameters; 2) gap fill the data, describe seasonal patterns of CO 2 exchange components and parameters, and estimate weekly and annual totals of gross primary production, total ecosystem respiration, and net ecosystem production (NEP); 3) compare CO 2 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.

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 BR 1 , 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 BR 2 ). 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 (BR 3 ) 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 BR 4 ) on a Zigweid fine-loamy soil ( (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 m 2 m −2 .
The BREB site in southeastern Wyoming (BR 5 ), 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). CO 2 exchange measurements at BR 5 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 EC 1 ) and moderately grazed (location EC 2 ) USDA Conservation Reserve Program (CRP) pastures at the Curtis Ranch near Briggsdale, Colorado, from 2004−2007, and on a pasture seeded to warm-season C 4 perennial bunch grass Caucasian bluestem (Bothriochloa bladhii [Retz] S. T. Blake) in the Texas High Plains (location EC 3 ) from 2010− 2011 (Rajan et al., 2013). Before 1987, the EC 1 (and EC 2 ) 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 wellestablished cover of the cool season C 3 grasses western wheat and needle-and-thread and the native warm season C 4 grasses buffalo grass, sideoats grama (Bouteloua curtipendula [Michx.] Torr.), and little Site EC 3 , 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). 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 PCPN h = 545 mm, sum of temperatures above 5°C T sum5 = 2 123-degree days) to the hottest and driest days in 2002 (PCPN h = 160 mm, T sum5 = 2 344-degree days). At the Wyoming site, lower temperatures (MAT = 6.7°C) and higher precipitation (PCPN h = 437 mm) in 1997 resulted in higher net production (NEP = 436 g CO 2 m −2 yr −1 ) than in 1998 (MAT = 7.8°C, PCPN h = 247 mm, NEP = 142 g CO 2 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 (PCPN h = 680 mm), followed by hot and dry conditions during the second (2011) year (MAT = 15.4°C, PCPN h = 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., 1997Falge 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 CO 2 fluxes (F) and the ancillary variables (incoming photosynthetically active radiation [Q], air temperature [T a ], 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 (P g ) and ecosystem respiration (R e ) components using the "lightsoil temperature-VPD" response method (Gilmanov et al., 2013(Gilmanov et al., , 2014Morgan et al., 2016). Following are the basic equations of the method: where α is the initial slope (apparent quantum yield), A max 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 k T are the coefficients of the exponential temperature response (r 0 = R e [0]), and the normalized VPD-response function φ(VPD,VPD cr ,σ VPD ) depends on two parameters: the critical value VPD cr , below which water deficit doesn't affect photosynthesis (φ= 1 for VPD ≤ VPD cr ), 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 r d was calculated as: where t 1 and t 2 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). Parameters of Eqs. [1-5] characterizing diurnal dynamics of the CO 2 exchange, such as apparent quantum yield α, photosynthetic capacity A max , 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 CO 2 m −2 s −1 ) of P g , R e , and F were calculated by Eqs. [1-5] using the diurnal data for meteorological drivers (Q, T a , T s , RH, VPD). The 24-hr integration of these rates provided the year-round daily (g CO 2 m −2 d −1 ) series of P g (t), R e (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][2][3][4][5] were also obtained. For days when identification of Eq. [1-2] parameters was not possible (mostly outside the growing season), the net CO 2 exchange for the day j was described as F(T s ) = − R e (T s ) with parameters of Eq.
[3] estimated from the subhourly (F, T s ) 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 Q a , ε phys = P g /Q a , and the ecological LUE coefficient, ε ecol = P g /Q, where Q is incident photosynthetically active radiation (Q a b Q). In contrast to ε phys characterizing mostly physiological abilities of plants to assimilate atmospheric CO 2 , ε ecol also incorporates ecologically significant information about structure and architecture of the plant canopy, as ε ecol = ε phys • f PARa (LAI), where f PARa (LAI) = Q a /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 lightuse 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 b 10 −6 ) between net primary production and ecological LUE. On these ecosystem comparisons, we are using the .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 T s ). ecological LUE coefficient for brevity denoted below as ε = P g /Q based on the incident PAR, Q.

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(Wylie et al., , 2007Gilmanov 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 BREB 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.

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 CO 2 exchange. Two major environmental response patterns of the CO 2 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 BR 2 site data for 2 June, 2003, which was a clear, sunny day (daily PAR total 51.9 mol m −2 d −1 , Q max = 2038 μmol m −2 s −1 ) with moderate temperatures (mean T a = 14.2°C, T a.max = 21.6°C) and low evaporative demand (VPD avg = 0.95 kPa, VPD max = 1.68 kPa). As shown in Figure 4A, on this day, CO 2 uptake of the SGS community is fairly well described by the univariate nonrectangular light-response function: where r d is average daily ecosystem respiration rate and other parameters as described above in Eqs. [1][2][3][4][5]. This equation describes a significant part of the variance of the CO 2 exchange (R 2 = 0.95) achieving the mean squared error of SE = 0.063 mg CO 2 m −2 s −1 . There is little difference between the morning and afternoon branches of the lightresponse curve for day of the year (DOY) 153. Predominance of radiation as a major driver of CO 2 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 CO 2 m −2 s −1 , R 2 = 0.96). Visually, the low significance of the VPD as a factor controlling the CO 2 uptake in this case is reflected by the fact that the black dots in Fig. 4B 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 , Q max = 2071.9 μmol m −2 s −1 ) with high temperatures (mean T a = 20.1°C, T a.max = 35.6°C) and strong evaporative demand (VPD avg = 2.5 kPa, VPD max = 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 r night for the 9-d window centered on DOY 182 showed a decrease of respiration with negative exponential temperature response coefficient k T = −0.058 (°C) −1 (see Fig. 5C). This suggests that the decrease of the CO 2 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. This conclusion is strongly supported by results of the nonlinear regression, which showed a highly significant (P b 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 k T = − 0.041 (°C) − 1 of low significance (P value 0.076). VPD-response function [2] with parameters VPD cr = 1.0 kPa and σ VPD = 3.16 kPa indicates a nearly fivefold decrease of CO 2 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 CO 2 exchange on DOY 182, 2003, the black dots in Fig. 5B (full model [1]-[4] predictions) markedly deviate from the VPD-independent surface F VPDavg (Q, T s ) = F(Q, T s , VPD avg ) generated from Eq. [2], with the VPD fixed at its mean value for day 182, VPD avg = 2.5 kPa, and are close to the gray dots (measurement data).
While plots similar to Figs. 4-5 were available for most of the days with valid {(Q, T s , 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.   (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 CO 2 m −2 s −1 for daily photosynthetic capacity A max 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 CO 2 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) (Zhang et al., 2012). 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, (C 4 ) and subdominant Pascopyrum smithii (C 3 ) 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.
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 ≤ r d ≤ 0.5 mg CO 2 m −2 s −1 , is approximately one-third the range of A max (Fig. 6C−6D). Pulses of intense CO 2 evolution observed in certain years (2000 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 (α, r d ) than observed at SGS proper sites (see Table 2). The highest values of parameters were achieved at the seeded Caucasian bluestem (C 4 ) pasture near Lockney, where in the rain-abundant year 2010, α max = 51 mmol mol −1 , ε max = 28.7 mmol mol −1 , and A max = 2.10 mg CO 2 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 T a indicates mean annual temperature; PCPN, annual precipitation; PCPN h , hydrologic year precipitation (October previous year to September current year); P g.max , maximum daily gross photosynthesis; P g.max.wk , maximum mean weekly gross photosynthesis; R e.max , maximum daily ecosystem respiration; R e.max.wk , maximum mean weekly ecosystem respiration; F max , maximum daily net CO 2 exchange; F max.wk , maximum mean weekly net CO 2 exchange; GPP, gross primary production; RE, ecosystem respiration; NEP, net ecosystem production. from the tallgrass ellipse, but it overlaps with the mixed-grass ellipse.
[1]-[5] are ecosystemscale (= canopy-scale as opposed to leaf-scale) parameters, meaning that, e.g., α is calculated per unit of incoming (not absorbed) PAR, and A max measures photosynthetic capacity per m 2 of ground surface (not per m 2 of leaf area, as in physiological photosynthetic capacity A max,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 A max 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 A max is proportional (directly proportional in monocultures) to the leaf area index LAI and A max,L , resulting in lower photosynthetic capacity of SGS canopies with LAI ≤ 1.5 m 2 m − 2 compared with tallgrass canopies with LAI ≥ 3 m 2 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 A max (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 P g (t), respiration R e (t), and net CO 2 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 CO 2 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 CO 2 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 (PCPN yr = 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 CO 2 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. Table 4 and Figure 9A show summary characteristics of the CO 2 exchange in ecosystems of the SGS ecoregion. Annual photosynthetic uptake varied greatly from b 1 000 to nearly 4 000 g CO 2 m −2 yr −1 , driven mostly by moisture availability. At the same time, respiration losses varied from 900 to 3 700 g CO 2 m −2 yr −1 , resulting in performance variation from sinks with NEP up to 700 g CO 2 m −2 yr −1 , to sources with NEP as low as −900 g CO 2 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 CO 2 -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 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).
(two-sided P values N 0.05) between means of gross primary production, ecosystem respiration, net ecosystem production, and ecophysiological parameters A max , r d.max , α max , and ε max (Fig. 10). The result remains unchanged (P values for differences of the means N 0.05) after exclusion of the outlier (site EC 3 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 b 0.05). This result is in agreement with the observation by Milchunas et al. (2008) that 50-yrlong 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 CO 2 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 BR 3 and BR 4 sites and the ungrazed and moderately grazed EC 1 and EC 2 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 P g for BREB sites are consistently higher than the curves for EC sites, while the respiration curve R e 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 CO 2 m −2 ) are higher than for the EC sites (3 676 and 3 982 g CO 2 m − 2 ), while the 3-yr RE cumulates for the BREB sites (3 792 and 3 789 g CO 2 m − 2 ) are lower than for the EC sites (4 970 and 4 920 g CO 2 m − 2 ) ( Table 5). This might be considered as Figure 10. Box-whisker plots of the CO 2 exchange components (A), maximum parameters of photosynthetic capacity and respiration rate (B), and quantum yield and light-use efficiency (C) of the ungrazed and grazed SGS ecosystems. 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 (BR 3 ), moderate continuously grazed (BR 4 ), ungrazed (EC 1 ), and moderately grazed (EC 2 ) 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. confirmation observations by Alfieri et al. (2009) that the BREB method generates higher net CO 2 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 b 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).
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.
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 P g (t) and R e (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 CO 2 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 (BR 3 and BR 4 ), 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 CO 2 exchange in rangelands of the SGS ecoregion, this study demonstrated considerable variability of daily maxima for photosynthesis 13 b P g b 41 g CO 2 m − 2 d − 1 , for respiration 13 b R e b 35 g CO 2 m − 2 d − 1 , and for net daily CO 2 flux 9 b F b 24 g CO 2 m − 2 d −1 . The annual ranges were 1 100 b GPP b 2 700 g CO 2 m − 2 yr −1 , 900 b RE b 3 000 g CO 2 m −2 yr −1 , and -900 b NEP b 700 g CO 2 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 CO 2 exchange were observed and the ecosystem response varied from a net sink of +700 g CO 2 m −2 yr −1 in 1999 with 487 mm precipitation to a net source of -900 g CO 2 m − 2 yr −1 in 2010 with 173-mm precipitation. Parameterization of CO 2 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 A max , and ecosystem respiration rate r d , including determination of the ranges of variation (19 b α b 51 mmol mol −1 , 7 b ε b 29 mmol mol −1 , 0.48 b A max b 2.1 mg CO 2 m −2 s −1 , 0.15 b r d b 0.49 mg CO 2 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 CO 2 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 P g and R e remaining within a ± 1.5 SD band around 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 GPP indicates gross primary productivity; RE, total ecosystem respiration; NEP, net ecosystem production. 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 A max ; (B) maximum light-use efficiency ε versus maximum daytime respiration rate r d . Black 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.
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 CO 2 exchange data, especially extensive in rangeland ecosystems, in ongoing efforts of comparative analysis, synthesis, and upscaling of rangeland CO 2 exchange and productivity.

Implications
The light-soil temperature-VPD-based method provides consistent partitioning of the raw flux tower CO 2 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 CO 2 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 CO 2 depending on weather. Quantification of the rates and parameters of the SGS ecosystems presented in this paper strengthens the empirical basis for sB60patiotemporal modeling and upscaling as tools for forage production and carbon management of the SGS rangelands.