Carbon response of tundra ecosystems to advancing greenup and snowmelt in Alaska
Study sites: site description and climatic limitIn this study, we focused on seven flux tower sites located in Alaska, United States, including six AmeriFlux sites and one site of the Korea Polar Research Institute (KOPRI) (Fig. S1A, Table S1). Over the seven study sites, the annual mean temperature was between −10.09 and −0.55 °C and the annual total precipitation ranged from 287 to 540 during 2001–2018 based on the North American Regional Reanalysis (NARR41, 0.3-degree resolution every 3 h). The annual mean temperature increased from 2001 to 2018 at all sites, with rates between 0.5 and 2.2 °C per decade (p 0.05 at all sites). Our study sites were mostly dominated by wet sedges, grasses, moss, lichens, and dwarf shrubs. For example, the dominant plants at the US-Atq site (at a higher latitude) are herbaceous sedges (Carex aquatilis, Eriophorum russeolum, and Eriophorum angustifolium) and shrubs (Salix rotundifolia), with abundant mosses (Calliergon richardsonii and Cinclidium subrotundum) and lichens (Peltigera sp.)11. At the KOPRI site (at a lower latitude), mosses (Sphagnum magellanicum, Sphagnum angustifolium, and Sphagnum fuscum), lichens (Cladonia mitis, Cladonia crispata, and Cladonia stellaris), and tundra tussock cottongrass (Eriophorum vaginatum) are abundant39. The active layer thickness is between 0.33 and 1.0 m, according to field data and radar-based estimates.Climatic limits imposed by temperature, water, and radiation were quantified following Nemani et al.42 at each site during the GS between 2001 and 2018 (Fig. S2) using the NARR data. In this study, we defined the GS as from May to Oct., early GS as between May and Jun., peak GS as between Jul. and Aug., and late GS as between Sep. and Oct. For a temperature limit scalar, the monthly mean temperature from −5 to 5 °C was linearly scaled between 100% (i.e., no growth) and 0% (i.e., no reduction in growing days). The monthly ratio of precipitation to potential evapotranspiration (PET by the Priestley-Taylor method43), ranging between 0 and 0.75, was linearly scaled from 100 to 0% as a water limit scalar. A radiation limit scalar was estimated as a 0.5% reduction in growing days for every 1% increase in monthly cloudiness above the 10% threshold (monthly cloudiness (n) was estimated44 as (R={R}_{0}(1-0.75{n}^{3.4})), where R and R0 are the monthly mean incoming radiation and clear-sky radiation45, respectively).The carbon flux response to climatic variations at each site was further analyzed using a forward stepwise multiple regression analysis11 between the NEE and meteorological variables (temperature, PAR, and VPD) using tower data during the GS. Interaction terms among the variables are also included to consider the convolved effects of the variables (Eq. (1)).$${Y}_{{{NEE}}}= {beta }_{0}+{beta }_{1}{X}_{{{T}}}+{beta }_{2},{X}_{{{VPD}}}+{beta }_{3}{X}_{{{PAR}}}+{beta }_{4}{X}_{{{T}}}* {X}_{{{VPD}}}+ldots$$$$qquad;;; {beta }_{5}{X}_{{{T}}}* {X}_{{{PAR}}}+{beta }_{6}{X}_{{{VPD}}}* {X}_{{{PAR}}}+{beta }_{7}{X}_{{{T}}}* {X}_{{{VPD}}}* {X}_{{{PAR}}}$$
(1)
where YNEE is the daily average NEE (µmol m−2 s−1), and XT, XVPD, and XPAR are daily average air temperature (°C), VPD (ha), and PAR (µmol Photon m−2 s−1), respectively. Regression coefficients (({beta }_{0},ldots ,{beta }_{7})), standard errors, significance (P-value), and R2 value of the final regression model are summarized in Table S2.MODIS: long-term trends of snowmelt and greenup timingsWe collected the gridded MODIS snow cover (MOD10A1.V00646 at a 500-m resolution every day) and phenology (MCD12Q2.V00634 at a 500-m resolution yearly) from the NASA Earthdata (https://earthdata.nasa.gov/). We estimated the snowmelt and snowpack timings at each site as the date when a logistic fit to the MODIS snow cover (quality flags of good and best) passed 0.1 each year. We rejected those snowmelt timings when the gaps in the daily MODIS snow cover were longer than 2 weeks around the time of the snowmelt event. The greenup and dormancy timings with a quality flag of best were taken from the MODIS phenology. Based on the spatial representativeness assessment (see Supplementary Note), we decided to use the snowmelt timing and greenup timing within a 1 × 1 pixel window.The significance of the long-term trends in greenup and snowmelt timings at each site was determined by Spearman’s rho and Mann-Kendall tests (Fig. S3). We further estimated the 95% confidence intervals of the trends from 3000 timing sets generated by bootstrap resampling from a normal distribution47 (mean equal to each greenup or snowmelt timing with three standard deviation set to 10 or 6.6 days, respectively, i.e., the root-mean-squared values between the ground data-based estimates and MODIS values in a 1 × 1 pixel window, Figs. S9 and S10).SGSI: snowmelt-growing season indexGrowing season index (GSI)48 is one of the novel phenology models49 and has been widely applied for the phenological representations of various ecosystems50,51. GSI is a product of three indices of climatic variables (Eq. (2), Fig. S5): daylength, VPD, and growing-degree-days (GDD)52. As a phenological measure for a given meteorological condition, we calculated the daily GSI for spring (from Jan. 1 to Jul. 31) and fall (from Aug. 1 to Dec. 31), respectively. For the spring-GSI, GDD is the degree sum when the daily mean temperature rises above −5 °C after Jan. 1. For the fall-GSI, GDD is the degree sum when the daily mean temperature falls below 20 °C after Aug. 1. We then revised the GSI by multiplying it by a snowmelt index (iS) and referred to this as the snowmelt-GSI (SGSI, Eq. (3), Fig. S5). This guarantees that vegetation greenup does not start unless snow is melted, even if the meteorological conditions are sufficient to trigger leaf-out. The iS was estimated to be 0 when the snow cover fraction (snowfac variable53 in ED2) was above 0.1 and 1 otherwise.$${{{GSI}}}={{iX}}_{1} times {{iX}}_{2} times {{iX}}_{3}$$
(2)
$${{{SGSI}}}={{{GSI}}} times {iS}$$
(3)
where iX (X1, X2, and X3 represent daylength, VPD, and GDD, respectively) is 0 (X ≤ Xmin), 1 (X ≥ Xmax), and (X − Xmin)/(Xmax − Xmin) otherwise. Xmax and Xmin are the maximum and minimum thresholds of each index, respectively. For the spring-GSI, Xmin was calculated as the minimum value among the values on the greenup day (from MCD12Q2.V006) for the study period of 2001–2018 at each study site, and Xmax was calculated as the minimum value among the values on the maturity day. For the fall-GSI, similarly, Xmin was the minimum value for the dormancy timings, and Xmax was the minimum value for the senescence timings. We incorporated GSI (or SGSI) into ED2 by multiplying it to the optimal value of leaf biomass on the day, where it operates as an upper limit of the leaf biomass.In this study, it was assumed that phenological stages are driven by meteorological conditions, not by other factors (e.g., no assumption of fixed phenological periods6,54,55). The development of a robust phenological model for the tundra ecosystem would be enabled by an increasing amount of ground-based phenology data (e.g., PhenoCam data37).Case study: flux data analysisThere were three sites where flux data is available for >5 years in Alaska; US-Atq site (flux data during 2004–2008 with delayed snowmelt in 2005), US-EML site (flux data during 2009–2017 with delayed snowmelt in 2017), and US-BZF (flux data during 2012–2018 with delayed snowmelt in 2017 and 2018). We first calculated two timings that are related to the meteorological conditions (0.1-GSI timing and half-max GSI timing51, Fig. S6) using the NARR data. The 0.1-GSI timing and half-max GSI timing were calculated on the day when the GSI passed 0.1 and the half-max value (i.e., 0.5), respectively, each year. To calculate two timings regarding the flux seasonal profile51 (source-sink transition timing and half-max productivity timing, Fig. S6), we used 30-min gap-filled FLUXNET201556 data (NEE and GPP; quality flags of measured or good) at the US-Atq site to calculate daily NEP (i.e., negative NEE) and daily GPP. At the US-EML and the US-BZF sites, we applied an open-source code called ONEFlux (Open Network-Enabled Flux processing pipeline for eddy-covariance data)57 using the ERA5 data (European Centre for Medium-Range Weather Forecast Reanalysis v558) which was downscaled with a quantile mapping method59. Using the gap-filled 30-min NEP and GPP data from the ONEFlux, we calculated the corresponding daily values and fitted a smoothing spline to the daily NEP and the daily GPP each year. The source-sink transition timing was defined as the day when the smoothing spline of the daily NEP passed zero in each year. The half-max productivity timing was set to the day when the smoothing spline of the daily GPP passed the half-max value in that year51.Further, we investigated whether the delayed snowmelt altered the relationships between meteorological conditions and the flux-threshold timings at each site based on (1) the correlation between the 0.1-GSI timing and the source-sink transition timing and (2) the correlation between the half-max GSI timing and the half-max productivity timing.ED2: model implementationWe used NARR data41 (0.3-degree resolution every 3 h; temperature, precipitation rate, pressure, v- and u-wind speed, downward longwave and shortwave radiation flux, and relative and specific humidity) for single-point ED2 implementation at each study site from 2001 to 2018. Vegetation structure (LAI, leaf and structural biomass, diameter at breast height, and population density) was initialized for each site by using the maximum annual LAI of cold-adapted shrubs and Arctic C3 grass from the Ent Global Vegetation Structure Dataset v1.0b (Ent-GVSD v1.0b) with the allometric equations in ED2. Ent-GVSD v1.0b provides plant functional types (from the MODIS land cover, MCD12C1.V00560) and maximum annual LAI values (from the MODIS LAI, MOD15A2.V00461,62) in subgrid cover fractions. We did not use canopy heights from Ent-GVSD v1.0b because of the absence of trees at the study sites. Soil texture (the ratio of sand:silt:clay) was set following the Harmonized World Soil Database v1.163 of the Food and Agriculture Organization of the United Nations (UN FAO). Soil carbon was initialized using the UN FAO Global Soil Organic Carbon Map64, and soil nitrogen was estimated using the soil C/N ratio of moist tundra (mean: 18.4)65.The prior distribution of each key variable was based on previous studies (Table S4), and 10,000 parameter sets were randomly generated from the prior distributions (the so-called Monte Carlo method). The best parameter set was selected based on statistical measures (r2 and root-mean-squared error) when compared to the data at the US-Atq site, i.e., NEP flux data for 2004–2006 and MODIS LAI data for 2003–2010 (MCD15A3H.V00666 at a 500-m resolution every 4 days) (Table S5). We then validated the performance of ED2 with this best parameter set by focusing on key ecosystem processes, such as NEP, ecosystem respiration, soil temperature, snowmelt timing, greenup timing, and the LAI at all sites (Table S5). The ED2 LAI was overestimated by 0.15–0.16 compared to the field-measured LAI values (Jul.–Aug. in 2006 at Barrow67 and Jun.–Aug. in 1996 at Toolik68).It is worth noting that the accuracy of the MODIS LAI has not been extensively evaluated at high latitudes because of limited ground measurements and few valid MODIS data points due to inadequate sun-sensor geometry, illumination conditions, and cloud contamination69,70. Furthermore, the heterogeneous landscapes of the region at the scale of remote sensing data (from hundreds of meters to a few kilometers) are also a major challenge that must be addressed before the data can be evaluated against ground measurements. According to the spatial representativeness assessment (see Supplementary Note, and Figs. S7 and S8), the landscapes around the flux towers generally have heterogeneity at a level similar to, or smaller than, the tower footprint size (200–300 m) during the early GS and peak GS in the MODIS 1 × 1 pixel window (i.e., 500 × 500 m2), but mostly higher than in the 3 × 3 pixel window (Table S6). This indicates that it is desirable to evaluate the MODIS 1 × 1 pixel LAI values against ground measurements, as both MODIS greenup and snowmelt timings agreed more with the ground data at the 1 × 1 pixel window scale than at the 3 × 3 pixel window scale (Figs. S9 and S10). A more thorough evaluation of both MODIS LAI data and ED2 LAI values is required in the coming years with the increase in ground data availability (e.g., National Ecological Observatory Network, NEON, LAI measurements).Correlation analysis: the effects of early or delayed snowmelt timingTo analyze the net and lagged effects of early or late snowmelt timing, it is necessary to constrain the contribution of interannual meteorological variation. Therefore, we compared only the years when meteorological conditions were similar, i.e., when the weekly mean GSI value was within one standard deviation of the weekly mean GSI during 2001–2018 (at the US-Hva site, weekly values during 1994–2018); meteorological conditions appeared similar for 8 or 9 years at each study site, except the US-BZF site, where the similarities were found for 10 years. We also limited the effect of greenup timing changes by excluding the years when greenup was earlier or later by one standard deviation of the mean of greenup timings during the study period. For the years satisfying the constraints, a least-squares linear regression was applied between the snowmelt timing change (deviation in snowmelt timing each year from the mean snowmelt timing) and the seasonal deviation (the difference of the seasonal mean from the mean value over the years) of each process from the ED2 results.To analyze the net and lagged effects of delayed snowmelt, we implemented the ED2 model in two schemes, (1) following the meteorologically-determined phenological index (i.e., the GSI, Eqs. (2) and (2) constraining leaf-out by snowmelt (i.e., the SGSI, Eq. (3)). For the years when unmelted snow delayed greenup, we took the difference between the modeled results (i.e., GSI and SGSI) for each process and applied a least-squares linear regression between the difference of each process and the delayed snowmelt days. More
