Impact of Shifts in Vegetation Phenology on the Carbon Balance of a Semiarid Sagebrush Ecosystem

Abstract

Dryland ecosystems are critical in regulating the interannual variability of the global terrestrial carbon cycle. The responses of such ecosystems to weather and environmental conditions remain important factors that limit the accurate projections of carbon balance under future climate change. Here, we investigated how shifts in vegetation phenology resulting from changes in weather and environmental conditions influenced ecosystem carbon cycling in one semiarid ecosystem in the Hanford area of central Washington, United States. We examined two years of measurements of the phenology camera, eddy covariance, and soil chamber from an upland semiarid sagebrush ecosystem. Both years had contrasting diel and seasonal patterns of CO${}_2$ fluxes, primarily driven by differences in vegetation phenology. The net ecosystem exchange of CO${}_2$ (NEE) and evapotranspiration (ET) in 2019 were enlarged by shifted vegetation phenology, as a cold and snow-covered winter and warm and dry winter in 2020 resulted in constrained magnitudes of NEE and ET during the summer months. The annual gross primary productivity (GPP) was much higher in 2019 than in 2020 (−211 vs. −112 gC m${}^{-2}$), whereas ecosystem respiration was comparable in these two years (164 vs. 144 gC m${}^{-2}$). Thus, the annual NEE in 2019 was negative (−47 gC m${}^{-2}$) with the sagebrush ecosystem functioning as a carbon sink, while the positive annual NEE in 2020 indicated that the sagebrush ecosystem functioned as a carbon source. Our results demonstrate that winter snowpack can be a critical driver of annual carbon uptake in semiarid sagebrush ecosystems.

1. Introduction

Dryland ecosystems constitute approximately 41% of the Earth’s land surface; thus, they have remarkable impacts on global carbon and water cycling [1,2]. Recent studies have illustrated that dryland ecosystem productivity plays an important role in regulating the trends and interannual variabilities of global terrestrial carbon sequestration [2,3,4]. These ecosystems, mainly including grasslands, shrublands, and savannas, are quite fragile and rather sensitive to changes in weather and climate conditions [5,6]. Enhanced warming and severe droughts have the potential to significantly alter the capacity of dryland ecosystems to sequester carbon [6,7]. Understanding how these ecosystems respond to weather and environmental variability is critical for improving our projections of global carbon sink–source dynamics.

Site studies and models have shown that water availability, indicated by soil moisture and precipitation, as well as evapotranspiration (ET), plays an important role in maintaining the productivity of dryland ecosystems (regarding changes that can significantly influence the structure and functioning of dryland ecosystems) [8,9,10]. In dryland regions, the annual net ecosystem exchange of CO${}_2$ (NEE) often switches between negative (carbon sink) and positive (carbon source) from year to year, mainly due to annual total precipitation crossing a threshold, which varies between different dryland ecosystems [9]. Vegetation structures have also been shown to influence soil water distribution [11,12]. Shrub species with rooting depths of several meters can vertically transport water from wetter to drier soil layers and maintain function during dry months [12]. Hence, woody species, such as shrubs in dryland ecosystems, usually respond slower to changes in environmental drivers than non-woody species, such as grasses, which rely heavily on shallow water availability [9].

The seasonal snowpack stores winter precipitation, which would replenish soil water storage by snow melting in the spring [13]. The snowpack also reduces soil water loss due to decreased evapotranspiration, resulting from an increase in the reflection of solar radiation by snow cover [14]. In the inland Pacific Northwest (iPNW), seasonal precipitation is the primary source of soil water availability for natural ecosystems with limited groundwater access [15], approximately 37% of which fall as snow [13]. However, the snowpack has shown strong interannual variability [16]. Given the widespread declining trends in both snow depth and snow cover duration in many regions of the world [17,18], examining the potential impacts of snowpack variability on dryland ecosystems is urgently needed for accurate projections of global carbon balance under future climate change.

In addition to the variabilities in NEE and gross primary production (GPP), variations in ecosystem respiration (R${}_{\mathrm{eco}}$) have also been associated with changes in soil water availability [19]. The temperature sensitivity of respiration in drylands is strongly controlled by soil moisture [19]. Dryland ecosystems commonly experience large pulses in R${}_{\mathrm{eco}}$ in response to heavy precipitation events, which account for a substantial portion of the ecosystem’s carbon balance [20,21]. In wet years, increases in GPP are often offset by a large R${}_{\mathrm{eco}}$ [9]. However, it is not yet clear how R${}_{\mathrm{eco}}$ responds to changes in soil water availability due to snow melt in semiarid ecosystems.

Phenology controls the seasonality of vegetation functioning and ecosystem feedback to changes in environmental drivers [22,23]. The increasing number of eddy covariance (EC) sites co-located with digital cameras (i.e., PhenoCam) for phenology monitoring have contributed to improving our understanding of the relationship between ecosystem productivity and phenology [24,25,26,27,28,29]. The green chromatic coordinate (G${}_{\mathrm{CC}}$) extracted from PhenoCam images is a plausible vegetation index used to track plant development throughout the season [25,26]. Past studies related to PhenoCam-derived G${}_{\mathrm{CC}}$ and EC-measured fluxes mainly focus on temperate and boreal forests [25], and only a few recent studies have focused on dryland ecosystems (e.g., [29]). In the iPNW, there remains a critical need to better understand the relationship between vegetation phenology and the carbon balance of dryland ecosystems.

The objective of this work was to evaluate the response of carbon fluxes to snowpack variability in a semiarid sagebrush ecosystem using a combination of eddy covariance, soil chamber, and phenology camera measurements during the years 2019 and 2020. In 2019, the semi-arid site experienced a cold and wet winter with a heavy snowpack that melted until the middle of March, while in 2020, the site experienced a warm and dry winter. As a result of the impact of snowpack in 2019, the seasonality of vegetation phenology was delayed, and we hypothesized that the ecosystem carbon balance (i.e., NEE, GPP, and R${}_{\mathrm{eco}}$) would also be altered. Specifically, our objectives were (1) to characterize how the snowpack influences vegetation phenology change and carbon uptake, as well as the diel and seasonal patterns of carbon fluxes in the semiarid ecosystem; (2) to compare the calculated R${}_{\mathrm{eco}}$ and measured soil respiration, and evaluate their performances in representing the respiration of dryland ecosystems.

2. Materials and Methods

2.1. Study Site

The experiment was conducted at one AmeriFlux site (US-Hn1; 46°24${}^{\prime }$32${}^{″}$N, 119°16${}^{\prime }$30${}^{″}$W; Figure 1) located in the Hanford area of central Washington, United States. Site information was documented in detail elsewhere [15,30]. Briefly, the area has a Mediterranean-type climate with generally cool, wet winters and hot, dry summers. The long-term mean annual precipitation is less than 200 mm, most of which occurs during winter and spring [15]. The site is characterized by a deep vadose zone in which the water input of the soil is primarily fed by precipitation. Vegetation is a mixture of scattered shrubs and short grasses (Figure 2). Shrub species include Artemisia tridentata (big sagebrush) and Chrysothamnus viscidiflorus (green rabbitbrush), with an average height of approximately 1.5 m. Grass species include Bromus tectorum (cheatgrass), Salsola kali (Russian thistle), Poa secunda (sandberg bluegrass), Pseudoroegneria spicata (bluebunch wheatgrass), and Stipa comate (needle-and-thread grass) [15]. The soil texture in the top 30 cm layer is loamy sand (5% clay, 11% silt, and 84% sand) with small rocks and gravel interspersed, and the soil texture in the layer from 30 to 45 cm is sand (2% clay, 4% silt, and 94% sand) [15,30]. An instrumented 10 m flux tower was installed at the site, allowing continuous eddy covariance and meteorological measurements [30]). The tower began to operate in January 2016 and a soil chamber system was installed in late 2018. For this study, we examined the measurements of the phenology camera, eddy covariance, and soil chamber systems collected during 2019 and 2020 with contrasting weather and vegetation conditions (Figure 2). The EC tower is located around 0.5 km from the river, at an elevation of 13 m higher than that of the river water surface [15]. Under the northwesterly and southeasterly dominated winds, the impacts of the river on the tower measurements are generally negligible, as indicated by the 90% contours of flux footprints (Figure 1). The flux footprint was determined using the tool developed by Kljun et al. [31].

2.2. Eddy Covariance Data

Eddy covariance (EC) measurements of energy, water vapor, and net ecosystem CO${}_2$ fluxes were made using a 3D sonic anemometer (Model CSAT3, Campbell Scientific, Inc., Logan, UT, USA) and an open-path gas analyzer (Model LI-7500A, LI-COR, Inc., Lincoln, NE, USA). The EC system was installed on the flux tower at 5 m above the ground, and the signals were sampled at a rate of 10 Hz and stored using a data logger (CR3000, Campbell Scientific, Inc.). Raw 10 Hz data were processed using EddyPro^® software (version 7.06, LI-COR, Inc.) to determine half-hourly fluxes of CO${}_2$ (NEE) and latent and sensible heat (LE and H). The 10 Hz data were filtered for out-of-range values and instrument failures and then despiked [32]. We applied the double rotation method [33] to the three-dimensional wind components. Then, block averaging was used to calculate the fluctuating parts for each 30-min interval. Covariances were calculated and corrected for the effects of high- and low-pass filtering [34,35] and air density fluctuations [36], respectively. According to Mauder and Foken [37], the corrected fluxes were quality-checked and flagged, and low-quality data were removed for this study. The “REddyProc” R package [38] was applied to determine the friction velocity (u${}_{*}$) threshold and fluxes under conditions of low turbulent mixing were also removed. These removed data were further gap-filled using a machine learning algorithm as described in Section 2.6.

2.3. Soil Chamber Data

Chamber measurements of soil respiration (R${}_{\mathrm{soil}}$) were made using an infrared gas analyzer (LI-8100A, LI-COR, Inc.) coupled to one automated long-term chamber (8100-104, LI-COR, Inc.; Figure 1). The soil chamber was placed over bare soil adjacent to patches of sagebrush and the soil collar was inserted into the soil with 2–3 cm exposed. The automated chamber was set to close for 3 min every 15 min to make measurements. CO${}_2$ concentration in the chamber was recorded every second over the sampling periods. R${}_{\mathrm{soil}}$ was calculated as an exponential fit of CO${}_2$ concentration versus time over each measurement interval using the SoilFluxPro™ software (version 4.0, LI-COR, Inc.). The average R${}_{\mathrm{soil}}$ of the two measurements within half-hour intervals was calculated to represent the 30-min mean soil respiration.

2.4. PhenoCam Data

To monitor vegetation phenology at the US-Hn1 site, we used a prescribed digital camera (NetCam SC IR, StarDot Technologies, Inc., Buena Park, CA, USA), configured and deployed according to a standard protocol of the PhenoCam network [25]. The digital camera was mounted on the tower at the same height as the EC system. Three-layer RGB images were recorded every 30 min and the “xROI” R package [26] was used to extract a time series of the green chromatic coordinate (G${}_{\mathrm{CC}}$) from the digital images across a delineated region of interest (ROI). The G${}_{\mathrm{CC}}$ is a commonly used metric of color information, which has been applied successfully in many ecosystems [28]:

$${\mathrm{G}}_{\mathrm{CC}}=\frac{{\mathrm{G}}_{\mathrm{DN}}}{{\mathrm{R}}_{\mathrm{DN}}+{\mathrm{G}}_{\mathrm{DN}}+{\mathrm{B}}_{\mathrm{DN}}},$$

(1)

where R${}_{\mathrm{DN}}$, G${}_{\mathrm{DN}}$, and B${}_{\mathrm{DN}}$ are the mean red, green, and blue digital numbers across the ROI, respectively. Digital images were screened to remove low-quality images, and daily G${}_{\mathrm{CC}}$ values were calculated from all available daytime images in a 3-day moving window [25,27]. In addition, we used the moderate-resolution imaging spectroradiometer (MODIS) enhanced vegetation index (EVI) to examine variations in vegetation activity in 2019 and 2020 compared to the long-term (2000–2018) average conditions. We obtained EVI data from the MOD13Q1 data product (https://lpdaac.usgs.gov/products/mod13q1v006/, accessed on 13 August 2021) at a 16-day temporal resolution and a 250-m spatial resolution at the tower location for the period 2000 to 2020.

2.5. Meteorological Data

A variety of meteorological instruments were also installed at the site. The meteorological data relevant to this study include solar radiation (R${}_{\mathrm{g}}$), air temperature (T${}_{\mathrm{air}}$), relative humidity (RH), wind speed (WS) and direction (WD), precipitation (P), soil temperature (T${}_{\mathrm{soil}}$) and volumetric water content (SWC) [15,30]. These data were sampled using the same data logger at 1 Hz and stored as 30-min averages. Furthermore, 15-min average meteorological data were obtained from two weather stations (Figure 1) located within 8 km of the EC tower; they were maintained by the Washington State University AgWeatherNet (https://weather.wsu.edu/, accessed on 13 August 2021). In the event of any gaps in the tower meteorological data, the 15-min data were averaged within half-hour intervals to fill the gaps [39]. Thus, the 30-min meteorological data used in this study were gap-free. Snow depth measurements were conducted at the Hanford Meteorological Station (HMS) Monitoring Network (http://www.hanford.gov/page.cfm/HMS, accessed on 7 April 2022). Furthermore, we obtained long-term data (1989–2020) at the location close to the US-Hn1 site from the PRISM Climate Group (http://prism.oregonstate.edu, accessed on 7 April 2022; spatial resolution of 4 km) to characterize weather and environmental conditions in 2019 and 2020 in the context of long-term climate conditions.

2.6. Flux Gap-Filling and Partitioning

In addition to the aforementioned quality assurance and control, field operations (e.g., instrument maintenance) and electrical and instrumental issues would also lead to data gaps. As a result, many EC sites lack approximately 20–60% of annual flux data [15,38,39,40]. At the US-Hn1 site, approximately 29%, 50%, 51%, and 40% of the total 30-min H, LE, NEE, and R${}_{\mathrm{soil}}$ data, respectively, were either removed due to low data quality or missed due to electrical and instrumental issues. Identifying a robust approach to filling the gaps in flux data has been an active research topic [39,41,42]. Recently, machine learning algorithms (for example, artificial neural networks and random forests) have been applied to fill the gaps in carbon fluxes [39,41,42]. While their performances largely degrade as the lengths of data gaps increase, the random forest (RF) algorithm has been found to outperform other approaches [39,41]. Here, we applied the RF algorithm to fill the data gaps in H, LE, NEE, and R${}_{\mathrm{soil}}$.

The RF algorithm [43] uses bootstrap aggregation and feature randomness to generate each individual tree, while the prediction results are calculated as an ensemble of many independent decision trees. The input variables used to train the RF algorithms included the above meteorological variables, vegetation index (EVI), three variables of the decimal day of the year, and sine and cosine functions representing seasonal changes [42]. The 16-day EVI data were resampled to 30 min data using cubic spline interpolation. In this study, we used the “randomForest” R package [44] and created 500 regression trees for each flux variable (H, LE, NEE, and R${}_{\mathrm{soil}}$). Furthermore, RF allows for an estimation of the input variables, is used to identify key drivers, and interprets the results [41].

After gap-filling, the “REddyProc” R package was applied to perform flux partitioning into GPP and R${}_{\mathrm{eco}}$,

$$\mathrm{GPP}=\mathrm{NEE}-{\mathrm{R}}_{\mathrm{eco}},$$

(2)

where R${}_{\mathrm{eco}}$ is calculated using fitted nocturnal temperature response functions (Reichstein et al., 2005) that describe changes to R${}_{\mathrm{eco}}$ (i.e., the nighttime method).

3. Results

3.1. Weather and Environmental Conditions

We used long-term climatology (1989–2018) and vegetation phenology (2000–2018) from the PRISM climate group and MODIS data product, respectively, as a reference to characterize the precipitation, air temperature, and vegetation activity in 2019 and 2020 (Figure 3). Annual precipitation in 2019 (215 mm) was 19 mm higher than the climatological mean annual total precipitation (196 ± 48 mm), while 2020 (150 mm) was 46 mm lower than the average. The annual mean air temperature in 2019 (11.1 °C) was 0.7 °C lower than the long-term average (11.8 °C), while 2020 (12.8 °C) was 1.0 °C warmer than the average. Compared to the 30-year average climatological conditions, 2019 was characterized by a cold and wet winter (Figure 3a,b). February and March of 2019 were 7.0 °C and 5.2 °C colder than the long-term averages of 3.5 ± 2.3 °C and 7.8 ± 1.3 °C, respectively. January–April 2019 were wetter than their long-term averages, with unusually high precipitation in February (64 mm vs. the climatological mean of 18 ± 11 mm). Based on the records from the HMS Monitoring Network and PhenoCam images, the area around the tower was covered by snow from 6 February to 15 March 2019 with an accumulated snow depth of approximately 0.8 m. The year 2020 was characterized by a warm and dry winter. The monthly mean air temperatures in January and February were 2.9 °C and 1.9 °C warmer than their climatological averages, respectively. The monthly precipitations in February through April were lower than their climatological averages with only a few centimeters of accumulated snow.

The long-term (2000–2018) phenology (Figure 3c) illustrated a well-defined annual cycle of greening and senescence of the sagebrush ecosystem around the tower location. Generally, the growing season was from February to June, with peak vegetation activity in April. Compared to long-term averages, the years 2019 and 2020 represented contrasting characteristics of vegetation phenology. In 2019, the vegetation greenup date was delayed, with vegetation activity peaking in May, and the monthly mean EVI in 2019 remained well above its averages from May through October. However, in 2020, the vegetation greenup date was slightly earlier than its average, with the peak vegetation activity occurring in March and the growing season was shorter compared to 2019.

Meteorological and PhenoCam measurements collected from the flux tower confirmed that 2019 and 2020 experienced contrasting meteorological and phenological conditions (Figure 4), especially in winter and early spring. In February and March of 2019, the snowpack provided a thermally stable and humid layer for the sagebrush ecosystem (i.e., approximately 0 °C of T${}_{\mathrm{soil}}$ and 12% of SWC at 5 cm, respectively), even though T${}_{\mathrm{air}}$ was well below 0 °C. The G${}_{\mathrm{CC}}$ value increased dramatically after snow melting and reached its maximum value at the end of April when the SWC decreased to approximately 8%. The grasses then faded and died in May, while the shrubs remained active until September. In 2020, with higher VPD and no snow cover, soil water loss continued from January through April. The G${}_{\mathrm{CC}}$ increased slightly in February and March with no clear peaks. In addition, there were a couple of rainfall events in the summer of 2020, resulting in a relatively wetter summer compared to the long-term average.

3.2. Comparison of the Diel Patterns of CO₂ Fluxes

To investigate the influence of different vegetation activities on carbon fluxes, we first compared the monthly diel variations of CO${}_2$ fluxes in 2019 and 2020 (Figure 5). NEE and GPP illustrated similar patterns with clear diel cycles in the growing season (i.e., March to July 2019 and February to April 2019, respectively) and constrained diel variations in the non-growing season. At the beginning of the growing season (e.g., March and April 2019; February and March 2019), NEE exhibited a U-shaped diel variability with the peak CO${}_2$ uptake around 11:00–12:00 Pacific Standard Time (PST). In the late months of the growing season (e.g., May to July 2019; April 2019), the magnitude of NEE reached its maximum value in the morning, then decreased slightly and reached a plateau in the early afternoon. The different diel patterns of NEE in the early and later growing seasons suggest different responses of vegetation types (grasses and shrubs) to environmental drivers. Regarding the monthly diel variations of ecosystem respiration, the calculated R${}_{\mathrm{eco}}$ followed similar diel patterns with minimum and maximum values appearing in the early morning and late afternoon, respectively. The difference between the maximum and minimum values was less than 0.3 µmol m${}^{-2}$ s${}^{-1}$. In months with a larger SWC and high temperature, R${}_{\mathrm{eco}}$ usually had larger values and larger amplitudes of diel variations.

3.3. Comparison of Seasonal Patterns of CO₂ Fluxes

To quantify the influence of different vegetation activities on the carbon balance of the ecosystem, we compared the seasonal patterns of CO${}_2$ fluxes in 2019 and 2020 (Figure 6). Both years showed distinct seasonal patterns of NEE. In 2019, the magnitude of monthly NEE was maximized in April (−28 gC m${}^{-2}$) and then gradually declined as the growing season progressed until July (Figure 6 and Table 1). The ecosystem was turned into a carbon source in August (13 gC m${}^{-2}$) and remained a weak source or sink in the following months. In 2020, the ecosystem carbon uptake was maximized in March (−24 gC m${}^{-2}$) and then declined dramatically in April (−8 gC m${}^{-2}$). The ecosystem functioned as a carbon source from May to December (Figure 6 and Table 1) due to lower vegetation activity and relatively higher ecosystem respiration compared to 2019.

The flux partitioning results illustrate that the different seasonal patterns in NEE in 2019 and 2020 were primarily driven by differences in GPP from February to October. In 2019, the magnitude of GPP was maximized in April (−57 gC m${}^{-2}$) and declined gradually in May through October to approximately −5 gC m${}^{-2}$. In 2020, the magnitude of the GPP was the largest in March (−34 gC m${}^{-2}$) and maintained relatively low values from May to October (Figure 6 and Table 1). The magnitude of total GPP was greater in 2019 than in 2020 (−211 gC m${}^{-2}$ vs. −112 gC m${}^{-2}$), whereas the annual total R${}_{\mathrm{eco}}$ of the two years was comparable (164 gC m${}^{-2}$ vs. 144 gC m${}^{-2}$). Overall, the ecosystem functioned as a carbon sink in 2019 and a carbon source in 2020 (−47 gC m${}^{-2}$ vs. 32 gC m${}^{-2}$).

3.4. Relationship between Vegetation Greenness and Ecosystem CO₂ Fluxes

Figure 7 shows the comparison of MODIS EVI and PhenoCam G${}_{\mathrm{CC}}$ and the distributions of NEE, GPP, and R${}_{\mathrm{eco}}$ with G${}_{\mathrm{CC}}$, respectively. The good fit between MODIS EVI and PhenoCam G${}_{\mathrm{CC}}$ (Figure 7) (R${}^2$ = 0.44; p < 0.001) suggests that PhenoCam can be a useful tool to continuously monitor vegetation phenology in dryland ecosystems. Note that the G${}_{\mathrm{CC}}$ data were screened to match the temporal resolution of the MODIS EVI data. In both years, G${}_{\mathrm{CC}}$ ranged from 0.325 to 0.375, with approximately 73% of the data less than 0.335. We fitted a third-degree polynomial function between NEE, GPP, R${}_{\mathrm{eco}}$, and G${}_{\mathrm{CC}}$, respectively, with the confidence level exceeding 99%, suggesting that G${}_{\mathrm{CC}}$ can be used as a metric to constrain NEE, GPP, and R${}_{\mathrm{eco}}$ at the dryland sagebrush ecosystem (Figure 7b–d). These results indicate that changes in vegetation phenology can significantly influence the carbon balance of the ecosystems in drylands. Long-term and accurate monitoring of vegetation phenology cannot only provide a ground-based reference for validating satellite measurements, but also provide a high-temporal resolution parameter for gap-filling the flux data. Therefore, complementing a low-cost PhenoCam at eddy covariance sites would greatly improve the quality of carbon flux data, especially for dryland ecosystems with sparse vegetation. Furthermore, with PhenoCam images, vegetation indices can be extracted and calculated for different vegetation types, which would be very helpful to assess the influence of different vegetation types on ecosystem CO${}_2$ fluxes in drylands.

3.5. Comparison of Soil and Ecosystem Respiration

As mentioned above, the small difference in the annual total ecosystem respiration despite distinct characteristics of meteorological conditions between both years inspired us to investigate the calculated R${}_{\mathrm{eco}}$ using the measured R${}_{\mathrm{soil}}$ as a reference. As shown in Figure 5d, the measured R${}_{\mathrm{soil}}$ also exhibited unimodal-shaped diel variations similar to R${}_{\mathrm{eco}}$ but with differences in the peak values. R${}_{\mathrm{soil}}$ increased quickly in the morning, reached its peak values at around noon, and then decreased gradually in the afternoon. The difference in the maximum ecosystem and soil respiration time indicates that there is still a need to improve the algorithms used for calculating ecosystem respiration, which is beyond the scope of the current study.

We also compared the daily sums of R${}_{\mathrm{eco}}$ and R${}_{\mathrm{soil}}$ over two years (Figure 8). We divided the data into four groups according to the ranges of daily mean SWC. Each group consists of daily sums/averages with SWC located within a range of 0–3, 3–6, 6–10, and 10–15%, accounting for 26.6, 41.4, 14.5, and 17.5% of daily data, respectively. For different SWC groups, we fitted a linear regression between R${}_{\mathrm{eco}}$ and R${}_{\mathrm{soil}}$. Under very dry conditions (0–3%), R${}_{\mathrm{eco}}$ was generally larger than R${}_{\mathrm{soil}}$ by approximately 43%, in agreement with the theoretical hypothesis since R${}_{\mathrm{soil}}$ is a component of R${}_{\mathrm{eco}}$. For soil moisture conditions within the range of 3–6%, R${}_{\mathrm{eco}}$ was quite close to R${}_{\mathrm{soil}}$, whereas, under relatively wetter conditions (i.e., SWC within 6–10% and 10–15%), R${}_{\mathrm{eco}}$ was in general smaller than R${}_{\mathrm{soil}}$, suggesting that under these conditions, the calculated R${}_{\mathrm{eco}}$ is generally underestimated the true ecosystem respiration. Furthermore, the seasonal patterns of R${}_{\mathrm{eco}}$ were quite similar to those of R${}_{\mathrm{soil}}$ in both years, and the monthly totals of R${}_{\mathrm{soil}}$ were generally lower than those of R${}_{\mathrm{eco}}$, except for several months, e.g., May 2019 (Figure 6 and Table 1). In general, the annual total soil respiration was 154 and 127 gC m${}^{-2}$ in 2019 and 2020, respectively, and the difference in the annual total R${}_{\mathrm{soil}}$ between both years was comparable to the difference in the annual total R${}_{\mathrm{eco}}$.

Table 1. Monthly average values of sensible heat flux (H, W m${}^{-2}$), evapotranspiration (ET, mm), net ecosystem exchange of CO${}_2$ (NEE, gC m${}^{-2}$), gross primary production (GPP, gC m${}^{-2}$), ecosystem respiration (R${}_{\mathrm{eco}}$, gC m${}^{-2}$), and soil respiration (R${}_{\mathrm{soil}}$, gC m${}^{-2}$) at the study site in 2019 and 2020, respectively.

	H		ET		NEE		GPP		R${}_{\mathbf{eco}}$		R${}_{\mathbf{soil}}$
	2019	2020	2019	2020	2019	2020	2019	2020	2019	2020	2019	2020
January	−1.0	−4.0	6.1	10.4	1.1	2.4	−11.5	−11.8	12.6	14.2	10.6	12.6
February	3.4	17.0	5.0	12.3	−3.1	−9.2	−7.2	−21.5	4.1	12.3	7.2	10.5
March	23.1	39.1	16.5	23.7	−6.7	−23.8	−17.4	−33.7	10.7	9.9	11.7	12.0
April	44.0	75.6	47.3	17.3	−28.3	−7.6	−57.2	−18.5	28.9	10.9	20.6	11.2
May	88.6	84.5	31.6	22.1	−17.6	16.1	−32.3	−3.5	14.7	19.6	17.7	14.9
June	92.1	97.1	20.7	23.2	−9.3	13.8	−23.8	−4.1	14.5	17.9	11.7	14.8
July	93.7	113.6	16.8	10.9	−5.1	13.3	−19.8	−0.9	14.7	14.1	11.9	9.0
August	78.4	89.8	17.3	7.9	12.7	7.4	−12.8	−2.6	25.4	10.0	20.3	7.1
September	50.1	51.3	13.2	7.9	4.0	3.2	−9.6	−4.8	13.6	8.0	12.4	6.0
October	28.1	29.3	8.4	8.2	4.5	7.6	−4.8	−1.1	9.3	8.7	11.9	6.9
November	9.7	5.2	5.3	8.3	−3.0	7.6	−7.5	−2.2	4.5	9.8	8.4	14.4
December	0.5	−0.4	3.9	4.5	3.9	1.2	−6.9	−7.0	10.8	8.2	9.1	7.9
Annual	42.4	50.0	191.9	156.6	−47.0	32.0	−210.9	−111.5	163.9	143.5	153.6	127.3

Table 1. Monthly average values of sensible heat flux (H, W m${}^{-2}$), evapotranspiration (ET, mm), net ecosystem exchange of CO${}_2$ (NEE, gC m${}^{-2}$), gross primary production (GPP, gC m${}^{-2}$), ecosystem respiration (R${}_{\mathrm{eco}}$, gC m${}^{-2}$), and soil respiration (R${}_{\mathrm{soil}}$, gC m${}^{-2}$) at the study site in 2019 and 2020, respectively.

	H		ET		NEE		GPP		R${}_{\mathbf{eco}}$		R${}_{\mathbf{soil}}$
	2019	2020	2019	2020	2019	2020	2019	2020	2019	2020	2019	2020
January	−1.0	−4.0	6.1	10.4	1.1	2.4	−11.5	−11.8	12.6	14.2	10.6	12.6
February	3.4	17.0	5.0	12.3	−3.1	−9.2	−7.2	−21.5	4.1	12.3	7.2	10.5
March	23.1	39.1	16.5	23.7	−6.7	−23.8	−17.4	−33.7	10.7	9.9	11.7	12.0
April	44.0	75.6	47.3	17.3	−28.3	−7.6	−57.2	−18.5	28.9	10.9	20.6	11.2
May	88.6	84.5	31.6	22.1	−17.6	16.1	−32.3	−3.5	14.7	19.6	17.7	14.9
June	92.1	97.1	20.7	23.2	−9.3	13.8	−23.8	−4.1	14.5	17.9	11.7	14.8
July	93.7	113.6	16.8	10.9	−5.1	13.3	−19.8	−0.9	14.7	14.1	11.9	9.0
August	78.4	89.8	17.3	7.9	12.7	7.4	−12.8	−2.6	25.4	10.0	20.3	7.1
September	50.1	51.3	13.2	7.9	4.0	3.2	−9.6	−4.8	13.6	8.0	12.4	6.0
October	28.1	29.3	8.4	8.2	4.5	7.6	−4.8	−1.1	9.3	8.7	11.9	6.9
November	9.7	5.2	5.3	8.3	−3.0	7.6	−7.5	−2.2	4.5	9.8	8.4	14.4
December	0.5	−0.4	3.9	4.5	3.9	1.2	−6.9	−7.0	10.8	8.2	9.1	7.9
Annual	42.4	50.0	191.9	156.6	−47.0	32.0	−210.9	−111.5	163.9	143.5	153.6	127.3

3.6. Influence of Environmental Drivers on Respiration in the Semiarid Ecosystem

By combining soil chamber measurements and flux partitioning results, we evaluated the response of respiration to environmental drivers in the semi-arid sagebrush ecosystem. For different SWC groups, we fitted the relationship between respiration (R${}_{\mathrm{eco}}$ and R${}_{\mathrm{soil}}$) and temperature (T${}_{\mathrm{air}}$ and T${}_{\mathrm{soil}}$) using exponential regression [45]. Our results show that for all four groups, the respiration generally increased with increasing temperature, although the increasing trends of respiration with increasing temperatures were different for the groups (Figure 9). The larger the SWC ranges, the faster the increasing trends (

Table 2. Coefficients (with 95% confidence bounds) of the fitted exponential relationships between respiration (R${}_{\mathrm{eco}}$ and R${}_{\mathrm{soil}}$) and temperature (T${}_{\mathrm{air}}$ and T${}_{\mathrm{soil}}$) in the form of y = Ae${}^{Bx}$.

		T${}_{\mathbf{air}}$		T${}_{\mathbf{soil}}$
		A	B	A	B
R${}_{\mathrm{eco}}$	0–3%	0.19 (0.11, 0.27)	0.03 (0.02, 0.05)	0.11 (0.06, 0.17)	0.04 (0.02, 0.05)
	3–6%	0.22 (0.18, 0.26)	0.04 (0.03, 0.05)	0.22 (0.18, 0.26)	0.03 (0.02, 0.04)
	6–10%	0.21 (0.17, 0.25)	0.07 (0.06, 0.08)	0.20 (0.16, 0.24)	0.05 (0.05, 0.06)
	10–15%	0.24 (0.21, 0.26)	0.09 (0.08, 0.10)	0.28 (0.24, 0.31)	0.06 (0.06, 0.07)
R${}_{\mathrm{soil}}$	0–3%	0.15 (0.09, 0.20)	0.03 (0.01, 0.04)	0.08 (0.05, 0.12)	0.04 (0.02, 0.05)
	3–6%	0.28 (0.21, 0.34)	0.02 (0.01, 0.03)	0.30 (0.23, 0.38)	0.01 (0.00, 0.02)
	6–10%	0.34 (0.27, 0.42)	0.04 (0.02, 0.05)	0.34 (0.26, 0.42)	0.03 (0.02, 0.04)
	10–15%	0.27 (0.22, 0.33)	0.08 (0.07, 0.10)	0.30 (0.25, 0.35)	0.06 (0.05, 0.07)

), suggesting that respiration is more sensitive to changes in temperature. Given the same temperature, respiration was higher for the larger SWC ranges, reflecting that soil water content affects ecosystem respiration. When soil or air temperature became lower, respiration from all groups converged to smaller values, though T${}_{\mathrm{air}}$ performed better than T${}_{\mathrm{soil}}$ in simulating respiration as T${}_{\mathrm{soil}}$ kept above 0 °C, leading to more data twinging around low T${}_{\mathrm{soil}}$ and thus a larger uncertainty. The change in the relationship between respiration and temperature for different SWC groups suggested that the calculated R${}_{\mathrm{eco}}$ reflected the impacts of soil moisture on respiration, although R${}_{\mathrm{eco}}$ was simulated as an exponential function of temperature in the nighttime method [46]. Therefore, we expect that at a daily scale, the calculated R${}_{\mathrm{eco}}$ was reliable to some extent in representing the influence of soil moisture and rain pulses even though the peaks in R${}_{\mathrm{eco}}$ and R${}_{\mathrm{soil}}$ at half-hourly scales were not matched to each other. However, due to the limited power supply provided by solar panels, only one soil chamber was installed for long-term measurements of soil respiration and, thus, it is hard to quantify the underestimation in respiration at this site.

4. Discussion

Here, we characterize and discuss the phenology of a semiarid sagebrush–grass ecosystem using MODIS EVI, a high-temporal resolution of G${}_{\mathrm{CC}}$ derived from PhenoCam, and measured carbon fluxes as well as the impacts of phenology shifts on the carbon balance. We observed contrasting seasonal variations in both PhenoCam-derived G${}_{\mathrm{CC}}$ and carbon fluxes in 2019 and 2020 (Figure 6 and Table 1). In general, seasonal variations in PhenoCam-derived G${}_{\mathrm{CC}}$ were in phase with those of carbon fluxes (Figure 7), consistent with the previous studies in other ecosystems (e.g., forests and grasses) [24,29]. We argue that the contrasting seasonal variations in G${}_{\mathrm{CC}}$ and carbon fluxes in 2019 and 2020 are primarily driven by the woody species at our study site [8]. The US-Hn1 site is characterized by relatively sparse shrubs (10%) and a large fraction of grasses [15]. Shrubs can utilize their vast rooting system to access groundwater deep within the soil [12,47], while grasses with dense shallow roots tend to use water in the topsoil layer [9]. In 2019, due to the presence of snowpack, the time of the G${}_{\mathrm{CC}}$ and carbon flux peaks were postponed, and the growing season for the shrubs was extended with available water from snowmelt, while in 2020, with warm and dry winter and early spring, the vegetation activity and carbon uptake peaked a couple of weeks earlier than in 2019 [48], and then decreased dramatically due to limited available water for both shrubs and grasses (Figure 4).

Environmental drivers play an important role in regulating seasonal variations of G${}_{\mathrm{CC}}$ and carbon fluxes (Figure 4 and Figure 6), though the role and importance of temperature and water availability vary across different phenological stages, especially for dryland ecosystems [8,9,29]. In winter and early spring, temperature is an important factor limiting ecosystem activities (e.g., photosynthesis and respiration) [48] (Figure 9). With a warm winter in 2020, the magnitudes of GPP and NEE in February were relatively larger than 2019 (Figure 6). In the late spring, summer, and autumn months, plant-available water is a limiting factor [8,9]. Due to the increase in solar radiation and air temperature as well as the associated increase in VPD, the increase in ET caused a rapid decrease in soil water availability in the topsoil layer (i.e., SWC) (Figure 4), and further resulted in the senescence of grasses and thus a reduction in carbon uptake [9], especially in 2020 (Figure 6). In 2019, we believe that plant-available water in a deep soil layer was replenished by snowmelt and, thus, shrubs maintained functioning during the dry months [12] (Figure 2). Note that the US-Hn1 site is an upland site with a deep vadose zone where the lateral groundwater–river water exchange is not a resource of plant-available water [49].

In summary, our results suggest that the presence of winter snowpack can be a critical driver for the annual carbon uptake of semiarid sagebrush ecosystems in snowy regions. A cold and wet winter that delays snow melt can result in shifts in the phenology of such ecosystems and further enhance the annual carbon uptake. However, global warming has been more prominent in winter than in other seasons, as indicated by the observed and modeled declining trends in both snow depth and snow cover duration in many regions worldwide [17,18,50]. How future climate change will affect the carbon balance in snowy semiarid regions still needs to be studied further due to the complex interactions among temperature, rainfall, snowpack, and evaporation [50].

5. Conclusions

In this study, we investigated the influence of vegetation phenology change on ecosystem carbon fluxes using two years of PhenoCam, eddy covariance, and soil chamber measurements from one semiarid ecosystem. Our results demonstrate that long, cold, and wet winters can be critical for the carbon balance in semiarid ecosystems. Due to low temperatures and the snowpack, vegetation activity was postponed in 2019, while snow-melted water was conserved in the deep soil layers and allowed the shrubs to maintain active until September. In contrast, the warm–dry winter in 2020 allowed the US-Hn1 site to retain a high NEE magnitude and ET in February and March when NEE was generally limited by low vegetation activity due to low temperatures. In dryland areas, vegetation activity is suspended by cold and wet winters, which can dramatically alter seasonal patterns of carbon balance.

Improving our understanding of the impacts of seasonal snowpacks on carbon cycling in semiarid ecosystems is crucial for better projection of how dryland ecosystems will respond to future changes in climate. The availability of water in such ecosystems is often a critical driver of carbon balance. The winter snowpack can fundamentally alter the temporal and spatial distributions of water availability and, thus, how dryland ecosystems will respond to changes in weather and environmental conditions. Long-term PhenoCam and eddy covariance measurements are required to comprehensively investigate the relationships between the semiarid ecosystem carbon balance and winter snowpack and climate conditions.