Standardized flux seasonality metrics: a companion dataset for FLUXNET annual product

Phenological events are integrative and sensitive indicators of ecosystem processes that respond to climate, water and nutrient availability, disturbance, and environmental change. The seasonality of ecosystem processes, including biogeochemical fluxes, can similarly be decomposed to identify key transition points and phase durations, which can be determined with high accuracy, and are specific to the processes of interest. As the seasonality of different processes differ, it can be argued that the interannual trends and responses to environmental forcings can be better described through the fluxes' own temporal characteristics than through correlation to traditional phenological events like bud break or leaf coloration. Here we present a global dataset of seasonality or phenological metrics calculated for gross primary productivity (GPP), ecosystem respiration (RE), latent heat (LE), and sensible heat (H) calculated for the FLUXNET2015 Dataset of about 200 sites and 1500 site years of data. The database includes metrics (i) on an absolute flux scale for comparisons with flux magnitudes and (ii) on a normalized scale for comparisons of change rates across different fluxes. Flux seasonality was characterized by fitting a single-pass double-logistic model to daily flux integrals, and the derivatives of the fitted time series were used to extract the phenological metrics marking key turning points, season lengths, and rates of change. Seasonal transition points could be determined with a 90 % confidence interval of 6–11 d for GPP, 8–14 d for RE, 10–15 d for LE, and 15–23 d for H. The phenology metrics derived from different partitioning methods diverged, at times significantly.

This Flux Seasonality Metrics Database (FSMD) can be accessed at the US Department of Energy's (DOE) Environmental Systems Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE, https://doi.org/10.15485/1602532; Yang and Noormets, 2020). We hope that it will facilitate new lines of research, including (1) validating and benchmarking ecosystem process models, (2) parameterizing satellite remote sensing phenology and PhenoCam products, (3) optimizing phenological models, and (4) generally expanding the toolset for interpreting ecosystems responses to changing climate.

1 Introduction

Phenology, the timing of life-cycle events and phases of plants and animals, and their relationship with the environment, especially climate (Lieth, 1974; Piao et al., 2019), is an important indicator of ecosystem dynamics. It is an integrating record of the effects of global warming and other environmental changes on biological processes (Noormets et al., 2009; Post and Stenseth, 1999; Weltzin et al., 2020). Current phenology studies focus primarily on structural changes such as bud break, flowering, leaf coloring, and leaf fall. However, the functional aspects of plant activities, although invisible, also provide quantitative measures of plant responses to changes in environmental conditions and underlie the structural changes (Fitzjarrald et al., 2001; Schwartz, 2003; Schwartz and Crawford, 2013). Ecosystem processes such as biogeochemical fluxes also show seasonal changes and can be decomposed to key transition dates and phase durations, which characterize the exchanges of mass and energy between plants and the environment and may exert mutual feedback (Baldocchi et al., 2001; Freedman et al., 2001).

Currently, the phenology datasets mainly have three sources: (i) ground-based observations of plant structural changes, (ii) camera-based observations of canopy reflectance (or near-surface remote sensing observations), and (iii) satellite-based observations of land surface reflectance. The ground-based phenology is a traditional but very useful method in phenology studies. For example, the USA National Phenology Network (USA-NPN) has collated observations of first bloom and first leaf of lilac and honeysuckle from the 1960s across the contiguous United States (Schwartz et al., 2012; CONUS, United States territory, not including Hawaii or Alaska; Betancourt et al., 2007; Glynn and Owen, 2015). The USA-NPN was established in part to assemble long-term phenology datasets for a broad array of species across the United States, which can be used to determine the extent to which species, populations, and communities are vulnerable to ongoing and projected future changes in climate (Glynn and Owen, 2015; Schwartz et al., 2012). The camera-based phenology observations such as the PhenoCam network (https://phenocam.sr.unh.edu/webcam/, last access: 21 March 2021; Richardson et al., 2018; Richardson, 2019) use high-resolution digital cameras to characterize canopy phenology through the color information from the images (Brown et al., 2016; Richardson et al., 2018). The remote sensing has been used to detect vegetation green-up and canopy development (Ganguly et al., 2010; Julien and Sobrino, 2009; Zhang et al., 2003, 2018). While the remote-sensing-based phenology product estimates transition dates from a continuous reflectance time series, it is compared against ground-based event data. The seasonality of ecosystem processes, including biogeochemical fluxes, can similarly be decomposed to identify key transition points and phase durations, which can be determined with high accuracy, and are specific to the processes of interest. As the seasonalities of different processes differ, it can be argued that the interannual trends and responses to environmental forcings can be better described through the fluxes' own temporal characteristics than through correlation to traditional phenological events like bud break or leaf coloration.

Therefore, the objective of this paper is to generate an objective and standardized flux seasonality metrics dataset, which can act as the companion dataset for the FLUXNET product. This study aims to develop a comprehensive framework for studying the seasonality of ecosystem processes systematically with eddy covariance flux measurements including gross primary productivity (GPP), ecosystem respiration (RE), latent heat (LE), and sensible heat (H). As ecosystem and Earth system models are increasingly tested and developed using the very high temporal and increasing spatial coverage of eddy covariance sites, the added information of explicit and standardized flux-specific transition times offers unprecedented opportunity to refine the process representations in models even further (Baldocchi et al., 2001; Falge et al., 2002; Noormets, 2009; Wofsy et al., 1993). The remainder of this paper is organized to the dataset description, the summary of estimating the seasonality metrics from idealized seasonal curves of different fluxes, description of the model performance, uncertainties of the flux seasonality metrics, and conclusions.

2 Data

FLUXNET is a global network of regional networks of eddy covariance sites that measure the exchange of CO₂, water vapor, and energy between vegetation and the atmosphere (Baldocchi et al., 2001; Baldocchi, 2008). Recently, harmonized data processing protocols have been developed (Pastorello et al., 2020), and the growing global coverage of these observations has become the de-facto ground-truthing tool for both mechanistic ecosystem models as well as global planetary circulation models (Baldocchi, 2003, 2020). The data include continuous (i.e., gap-filled) measurements of net ecosystem exchange of CO₂ (NEE), latent and sensible heat fluxes (LE and H), and microclimate data (air temperature, humidity, wind speed and direction, solar radiation, soil temperature, and soil water content), all at a 30 min time step. Estimates of canopy photosynthesis and ecosystem respiration, derived from the data using an empirical model, are also typically available. The data undergo quality assurance, and missing half-hourly averages are filled using standardized methods to provide continuous data records (Foken and Wichura, 1996; Mauder et al., 2008; Pastorello et al., 2019). The current study uses the FLUXNET2015 Dataset (https://fluxnet.fluxdata.org/, last access: 21 March 2021), which includes over 200 sites and around 1500 site years of data (Fig. 1). The gap-filled 30 min data series of fluxes and micrometeorological conditions were aggregated to daily totals. The example sites for each biome were selected based on the following boundary conditions: (1) distinct seasonality of all fluxes; (2) data coverage of observed and high-quality gap-filled data >75 % (defined by variable-specific data quality flags in the FLUXNET database; Reichstein et al., 2005). Therefore, the coverage of different fluxes is different, in which GPP has the highest coverage and H has the lowest coverage. The final dataset included 170 sites and 1035 site years for GPP, 173 sites and 1013 site years for RE, 130 sites and 773 site years for LE, and 110 sites and 446 site years for H. Even though the FLUXNET sites are mainly distributed in the Northern Hemisphere and temperate ecosystems, they still have high spatial and temporal representativeness (Yu et al., 2019).

Figure 1 Global distribution of eddy covariance flux sites included in this study, color coded by their International Geosphere-Biosphere Programme (IGBP) biome classes.

Although there is a broad agreement between different flux partitioning approaches (Moffat et al., 2007), and many approaches have converged recently, the current FLUXNET data product still includes a couple of alternative estimates of derived fluxes (RE, GPP). The latest interpretation of respiration fluxes, in particular, is that the night- and daytime estimates may represent different facets of the same phenomenon (Keenan et al., 2019). Different researchers may choose different partitioning approaches for different purposes. Hence, the phenological metrics dataset includes metrics for all biogeochemical fluxes reported in the FLUXNET dataset. Debating the strengths of different partitioning approaches, and indeed the quality of the underlying dataset, is beyond the scope of this study. We assume that the data reported has passed certain thresholds, even though some additional screening has been necessary, and not all site years are of sufficient coverage and quality to estimate the flux seasonality metrics. The metrics are reported both in absolute flux scale (to allow comparisons against commonly reported values) and in relative, normalized scale (to allow comparisons of development rates among different fluxes).

3 Method

3.1 Phenology metrics and uncertainties from flux observation

The seasonal dynamics of the ecosystem fluxes generally have five distinctive phases, which results from the interaction between the inherent biological and ecological processes and the changes in environmental conditions, and reflects the unique functioning of plant community at different stages of the growing season (Gu et al., 2009). The five phases are (1) pre-phase, baseline dormant season flux, before leaf development; (2) flux development period, a rapid increase in flux rate, concurrent or immediately following leaf emergence and expansion; (3) peak flux period, a relatively steady stage in the middle of the growing season; (4) flux recession period, a rapidly declining stage to the baseline; and (5) termination phase, the onset of a new dormant season, following leaf senescence and abscission. Sites with non-standard flux seasonality, where the R² of model fit was below 0.75, were filtered out during general model fit assessment. Most filtering related to poor quality gap filling with obviously distorted seasonality, but seven sites with multiple peaks during the same year were also excluded from the current study. Southern Hemisphere sites were analyzed by shifting the calendar year cutoffs by 180 d.

Here, we fit the single-pass double-logistic model as first described by Gu et al. (2009) to fit the flux time series. This function exhibits broad structural flexibility and robust convergence, both of which are important for automated processing. The temporal variation in eddy-flux data for an entire growing cycle can be modeled using the function

$$\begin{array}{}\text{(1)}& F_{\mathrm{m}}\left(t\right)=f_{\mathrm{0}}+{\displaystyle \frac{a_{\mathrm{1}}}{{\left(\mathrm{1}+e^{-b_{\mathrm{1}}\left(t-t_{\mathrm{1}}\right)}\right)}^{c_{\mathrm{1}}}}}-{\displaystyle \frac{a_{\mathrm{2}}}{{\left(\mathrm{1}+e^{-b_{\mathrm{2}}\left(t-t_{\mathrm{2}}\right)}\right)}^{c_{\mathrm{2}}}}},\end{array}$$

where F_m is the flux value in a given day of year (DOY) t, f₀ is the dormant season base flux, and a₁ and a₂ are parameters about the flux magnitude. Parameters b₁, b₂, c₁, and c₂ are related with the transitions or curvature parameters. The model was fit to daily integrated fluxes, following iterative procedures:

• a.

  Fit Eq. (1) to the flux time series and calculate a predicted value for each DOY.

• b.

  For each point in the time series, compute the ratio of the observed to predicted flux.

• c.

  Conduct the Grubbs test (Grubbs, 1969) to identify outliers in the obtained ratios.

• d.

  If an outlier is detected, remove this outlier and go to step c.

• e.

  If no outliers are found, remove the data points whose ratios are more than 1 standard deviation (1 s) below the mean ratio.

• f.

  Fit Eq. (1) to the time series of the daily flux measurements.

The DOYs at which the fitted logistic curve showed characteristic curvature changes were identified with the formula shown in Table 1 derived analytically from the seven parameters of Eq. (1) corresponding to the minimum and maximum values of the first and second derivatives. The first derivative of Eq. (1) is given by

$$\begin{array}{}\text{(2)}& F_{\mathrm{m}}^{\prime }\left(t\right)={\displaystyle \frac{a_{\mathrm{1}}{c}_{\mathrm{1}}{e}^{-\left(b_{\mathrm{1}}\left(t-t_{\mathrm{1}}\right)\right)}}{b_{\mathrm{1}}(\mathrm{1}+e^{-b_{\mathrm{1}}\left(t-t_{\mathrm{1}}\right)}{)}^{c_{\mathrm{1}}+\mathrm{1}}}}-{\displaystyle \frac{a_{\mathrm{2}}{c}_{\mathrm{2}}{e}^{-\left(b_{\mathrm{2}}\left(t-t_{\mathrm{2}}\right)\right)}}{b_{\mathrm{2}}(\mathrm{1}+e^{-b_{\mathrm{2}}\left(t-t_{\mathrm{2}}\right)}{)}^{c_{\mathrm{2}}+\mathrm{1}}}}.\end{array}$$

The first derivative can be considered as the rate of change of the flux. The maximum of the first derivative occurs early and the minimum late in the growing season (Fig. 2). The day on which the maximal growth rate of each flux occurs is termed the “midpoint of flux development period” (DOY_MFD; point B), the day on which the minimal growth rate occurs is termed the “midpoint of flux recession period” (DOY_MFR; point E), and the interval between these two transition dates is termed the “length of flux midpoint” ($L_{\text{FM}}=E-B$).

Figure 2 An example of the seasonal dynamics of gross primary productivity (GPP), and metrics of transition points of different phases derived from the extremes of the first ($F_{\mathrm{m}}^{\prime }\left(t\right)$) and second ($F_{\mathrm{m}}^{\prime \prime }\left(t\right)$) derivatives of the fitted logistic function (Eq. 1). For visual clarity, the scales of the first and second derivatives are enhanced 20-fold and 200-fold, respectively (orange and purple lines). The blue line indicates the double-logistic model (Eq. 1) fitted to the observed flux time series (black dots). The slope of the green dash lines indicates the rate of change during the flux development, peak flux, and flux recession period. The phenological transition points are marked with the vertical dashed lines, and the bootstrap estimates of 90 % confidence intervals of these metrics are indicated with the horizontal red error bars about the seven key transition points.

The second derivative of Eq. (1) is given by

$$\begin{array}{ll}\text{(3)}& {\displaystyle }{F}_{\mathrm{m}}^{\prime \prime }\left(t\right)=& {\displaystyle }{\displaystyle \frac{a_{\mathrm{2}}{c}_{\mathrm{2}}{e}^{-\left(b_{\mathrm{2}}\left(t-t_{\mathrm{2}}\right)\right)}}{b_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}+e^{-\left(b_{\mathrm{2}}\right(t-t_{\mathrm{2}})}{)}^{c_{\mathrm{2}}+\mathrm{1}}}}{\displaystyle }& {\displaystyle }-{\displaystyle \frac{a_{\mathrm{1}}{c}_{\mathrm{1}}{e}^{-\left(b_{\mathrm{1}}\left(t-t_{\mathrm{1}}\right)\right)}}{b_{\mathrm{1}}^{\mathrm{2}}{\left(\mathrm{1}+e^{-\left(b_{\mathrm{1}}\right(t-t_{\mathrm{1}})}\right)}^{c_{\mathrm{1}}+\mathrm{1}}}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{a_{\mathrm{1}}{c}_{\mathrm{1}}{e}^{-\left(\mathrm{2}{b}_{\mathrm{1}}\left(t-t_{\mathrm{1}}\right)(c_{\mathrm{1}}+\mathrm{1})\right)}}{b_{\mathrm{1}}^{\mathrm{2}}{\left(\mathrm{1}+e^{-\left(b_{\mathrm{1}}\right(t-t_{\mathrm{1}})}\right)}^{c_{\mathrm{1}}+\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{a_{\mathrm{2}}{c}_{\mathrm{2}}{e}^{-(\mathrm{2}{b}_{\mathrm{2}}\left(t-t_{\mathrm{2}}\right)\times (c_{\mathrm{2}}+\mathrm{1}\left)\right)}}{b_{\mathrm{2}}^{\mathrm{2}}(\mathrm{1}+e^{-\left(b_{\mathrm{2}}\right(t-t_{\mathrm{2}})}{)}^{c_{\mathrm{2}}+\mathrm{2}}}}.\end{array}$$

The second derivative can be considered as the rate of the growth rate of the flux. The spring and fall maxima of the second derivative mark the start of flux development period (DOY_SFD; point A) and end of the flux recession period (DOY_EFR; point F), whereas the minima mark the transition from the flux development to peak flux period ($\text{point C}={\text{DOY}}_{\text{EFD}}={\text{DOY}}_{\text{SPF}}$) and the transitions from the peak flux to flux recession period ($\text{point D}={\text{DOY}}_{\text{EPF}}={\text{DOY}}_{\text{SFR}}$). Periods between AC, CD, and DF mark the length of flux development, peak flux, and flux recession periods (L_FD, L_PF, and L_FR, respectively). We also calculated the peak flux (Fmax), date of peak flux (DOY_Fmax), the rate of the flux development period (R_FD;l₁), the rate of peak flux period (R_PF;l₂), and the rate of the flux recession period (R_FR;l₃). Period AF is the length of active season (L_AS).

The fitted daily fluxes observations were normalized to 0–1, and then the rates of change were calculated (namely NR_FD, NR_PF, and NR_FR), which can be used for comparison of development rates among different fluxes.

All the seasonality metrics were summarized in Table 1.

Table 1 Seasonality metrics estimated for biogeochemical fluxes including gross primary production, ecosystem respiration, latent heat, and sensible heat.

Note that “–” means units vary for phenology metrics of different fluxes.

3.2 Evaluation of parameter quality

3.2.1 Model fit statistics

The fit of Eq. (1) to flux time series was characterized through the coefficient of determination (R²), root mean square error (RMSE), empirical BIAS, and agreement index (AI).

The R² value of a regression is a measure of the portion of the variance of the dependent variable accounted for by the explanatory variables and characterizes the goodness of fit of the fitted model,

$$\begin{array}{}\text{(4)}& R^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n\left(F_{\mathrm{o}}\right(t)-\overline{F_{\mathrm{o}}}{)}^{\mathrm{2}}-{\sum }_{t=\mathrm{1}}^n(F_{\mathrm{o}}\left(t\right)-F_{\mathrm{m}}\left(t\right){)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^n\left(F_{\mathrm{o}}\right(t)-\overline{F_{\mathrm{o}}}{)}^{\mathrm{2}}}},\end{array}$$

where F_m and F_o are the predicted and observed values, respectively, $\overline{F_{\mathrm{o}}}$ is the mean value of the observations, and n is the sample size or the number of days in the year.

The RMSE was estimated as the square root of the mean value of the squared residuals:

$$\begin{array}{}\text{(5)}& \text{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n\left(F_{\mathrm{o}}\right(t)-F_{\mathrm{m}}(t){)}^{\mathrm{2}}}{n}}}.\end{array}$$

The BIAS was calculated as the mean value of the model's residuals:

$$\begin{array}{}\text{(6)}& \text{BIAS}={\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n\left(F_{\mathrm{o}}\right(t)-F_{\mathrm{m}}(t\left)\right)}{n}}.\end{array}$$

The agreement index (AI; Willmott, 2013) provides a measure of the relative error in model estimates, combining the information contained in the correlation coefficient (R) and RMSE, and is popular in model assessments and calibration (Gu et al., 2002; Zhou et al., 2016). It is calculated as

$$\begin{array}{}\text{(7)}& \text{AI}=\mathrm{100}-\mathrm{100}{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n\left(F_{\mathrm{o}}\right(t)-F_{\mathrm{m}}(t){)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^n{\left(\left|\left(F_{\mathrm{m}}\right(t)-\overline{F_{\mathrm{o}}}\right|+\left|F_{\mathrm{o}}\left(t\right)-\overline{F_{\mathrm{o}}}\right|\right)}^{\mathrm{2}}}}.\end{array}$$

AI is dimensionless and ranges from 0 (complete disagreement) to 100 (perfect fit). The AI is also sensitive to differences between observed and modeled means (Willmott, 2013). Thus, the AI is well suited for comparing model fits across different biomes and climates.

3.2.2 Uncertainties estimation

The uncertainties in the flux seasonality metrics estimates arise from two sources: (i) the day-to-day variability of fluxes, particularly during the transition periods, that affect the overall goodness of fit of Eq. (1) (see Sect. 3.2.1) and (ii) the consistency of change in climatological drivers during the transition periods that can manifest as early or late cold or warm spells providing conflicting signals to plant development and can affect specific metrics without affecting others. The overall model fit statistics can be used to identify the suitability of different data sources for flux seasonality assessment, but they are not good indicators of the confidence in specific seasonality metrics. Assessing the quality of the underlying flux data is beyond the scope of the current study, and all reported flux values are assumed “true” and the best possible estimates. The uncertainties in the seasonality metrics were estimated similar to Elmore et al. (2012), using Monte Carlo bootstrapping (Efron, 1979). Bootstrapping is a statistical procedure that resamples a single dataset to create many simulated samples. The advantage of the bootstrapping is that parameters can be estimated without assumptions about the normal distribution and using also a small sample size. The distribution of parameter estimates for these bootstrap models provides valuable information about parameter uncertainty and correlation that is free of assumptions about the underlying data distributions. In this study, random uniform sampling with replacement was conducted 500 times for each site year, and the seasonality metrics were estimated for each iteration of the bootstrapped dataset (Elmore et al., 2012; Klosterman et al., 2014). The 5th and 95th percentiles of the 500 bootstrapped phenology metrics estimates were taken as the confidence interval of the mean estimated from the original dataset. The uncertainties of the seven key transition dates are shown in Fig. 2, with red horizontal error bars to indicate the uncertainty intervals.

4 Results

4.1 Model fit statistics

The double-logistic model (Eq. 1) captured the temporal dynamics of widely divergent flux time series (Fig. 3). Although the model fit statistics do not directly translate to the quality of the seasonality metrics estimates (see Sect. 3.2.2), the general fit statistics deserve a brief review. Table 2 shows the overall performance of the fitted model for the different fluxes. The primary explanatory variable behind the fit statistics, as well as the differences between different fluxes, was the range of day-to-day variability in the flux time series. For example, H was generally more variable than LE, both of which are much more variable than RE and GPP, which resulted in lower fit statistics for H and LE than for GPP and RE (Table 2, Fig. 3). The model fit statistics reported in Table 2 were largely consistent in ranking the goodness of fit among biomes. Biomes with well expressed seasonal flux magnitude differences (mixed forest, deciduous broadleaf forest, evergreen needleleaf forest, croplands) exhibited consistently higher fit statistics than biomes with weak seasonality of fluxes (CSH and EBF). The latter also exhibited relatively greater day-to-day variability of fluxes, resulting in lower fit statistics. For GPP, the fit statistics were practically indistinguishable for time series partitioned based on daytime and nighttime partitioning methods (GPP-DT: AI=98.81, R²=0.96; GPP-NT: AI=98.73; R²=0.95), whereas for RE the nighttime partitioning method showed a marginally better fit than the daytime method (AI=98.66 and R²=0.95 versus AI=97.82 and R²=0.92).

Figure 3 Examples of the seasonal dynamics of different fluxes for 10 sites representative of different biomes. One biome, open shrubland, was left off because of space limitations on a single page. The blue line indicates the double-logistic model (Eq. 1) fitted to the observed flux time series (black dots). The phenological transition points are marked with the vertical dashed lines, and the bootstrap estimates of 90 % confidence intervals of these metrics are indicated with the horizontal red error bar for corresponding transition points. Note that the confidence intervals are not always symmetrical to the best estimate.

Table 2 Model fit statistics to different flux time series in different biomes.

R²: coefficient of determination; RMSE: root mean square error; AI: agreement index; GPP: gross primary production; RE: ecosystem respiration; H: sensible heat flux; LE: latent heat flux; DT: daytime partitioning method; NT: nighttime partitioning method; CORR: corrected H_F_MDS or LE_F_MDS; MDS: marginal distribution sampling; DBF: deciduous broadleaf forests; EBF: evergreen broadleaf forests; ENF: evergreen needle leaf forests; CRO: croplands; OSH: open shrublands; CSH: closed shrublands; GRA: grasslands; MF: mixed forests; SAV: savannas; WSA: woody savannas; WET: permanent wetlands.

4.2 Uncertainties of seasonality metrics

The uncertainties of individual flux seasonality metrics (Table 1) were estimated as the 5th and 95th percentiles of 500 Monte Carlo bootstrapping samples ranging from about a week to several weeks, and the uncertainties of phase durations tended to generally be greater than those of individual transition dates. Generally, the uncertainties were the lowest for the phenology metrics of GPP and highest for H (shown as horizontal red lines on Fig. 3). The average uncertainties of transition dates ranged from 6–11 d for GPP, 8 to 14 d for RE, 10 to 15 d for LE, and 15 to 23 d for H. The average uncertainties of duration length ranged from 12–20 d for GPP, 14 to 23 d for RE, 16 to 25 d for LE, and 23 to 32 d for H.

For all fluxes, deciduous broadleaf forest always showed the lowest uncertainties among all biomes. Uncertainties of flux development midpoints were always lower than those of start and end points. Meanwhile, the uncertainties of duration dates are larger than those of transition dates, indicating the compounding effect of the uncertainties of the start and end dates of the active season. The length of the dormant season also affects the uncertainties: the longer the dormant season, the lower the uncertainties.

4.3 Alternative data sources

4.3.1 Daily peak versus daily total fluxes

Even with standardized data sources like FLUXNET2015 Dataset, some user discretion of data aggregation remains, which may affect the reproducibility of different analyses. We mentioned earlier that the FLUXNET2015 Dataset contains GPP and RE estimates from different gap-filling algorithms. In addition, the daily time series can be assembled as either daily integrated fluxes (more common) or as daily peak fluxes, as proposed by Gu et al. (2009), reasoning that daily peaks may be less sensitive to cloudiness and thus better capture the seasonal dynamics of flux capacity. To test this line of argumentation, and identify the best representation of flux seasonality, we started the current study by comparing the model fit statistics for either of these data types. Except for H, all other fluxes had higher model fit statistics with daily integrated than daily peak fluxes (Fig. 4). Although only R² values are shown, all other fit statistics confirmed the same pattern (Table 2). Therefore, we adopted the daily integrated fluxes for the following phenology metrics dataset generation. Although the day- and nighttime partitioning methods yielded different flux estimates, and different model fits, choosing the best between them is not appropriate at this point. Both partitioning methods have their uses, and their respective merits have not been conclusively proven. The differences between the seasonality metrics of each dataset will be discussed in Sect. 4.2.3.

Figure 4 Comparison of the coefficients of determination (R²) of fitted seasonality curves between daily integrated fluxes (x axis) and daily peak fluxes (y axis) for GPP (a, b), RE (c, d), LE (e), and H (f). Both daytime (DT) and nighttime (NT) flux partitioning models are shown for GPP and RE.

Although daily flux totals were identified as the preferred scalar for deriving seasonality metrics from, we will report here the differences between the metrics estimated from the daily peak ecosystem respiration and daily integrated ecosystem respiration data (Fig. 5). This can be viewed as the minimum methodological uncertainty in a best-case scenario in the sense that the correlation between the two sets of metrics was much higher (average R²=0.81; Fig. 5) than between other sources of variability (e.g., the daytime and nighttime partitioning models resulted in seasonality metrics with an average R²=0.27; Fig. 7). The differences between the metric derived from daily integrated and peak flux values are generally smaller than the confidence intervals of individual estimates, except in WSA and SAV. Annual peak flux from daily peak ecosystem respiration exhibited consistently greater values than daily integrals (Fig. 5o), as would be expected.

Figure 5 The scatter plots of different phenology metrics from RE-NT daily integrated data and RE-NT daily peak data.

4.3.2 Comparison of different flux partitioning methods

The significance of the assumptions made by partitioning methods to fill flux time series has been emphasized from the perspective of flux integrals (Kruijt et al., 2004). Here we show that the choice of the partitioning model can also affect the seasonality of daily integrated fluxes and thus the seasonality metrics (Fig. 6). We report here the differences between the daytime and nighttime models of flux partitioning as exemplified by the RE and GPP time series, but the lessons apply for all partitioning approaches. Most importantly, mixing of time series filled with different models should not be done.

Figure 6 Seasonal dynamics of RE from daytime and nighttime partitioning methods for 11 sites representative of different biomes.

Figure 7 The scatter plots of seasonality metrics from RE data using daytime and nighttime partitioning methods.

Ecosystem respiration (RE)

The nighttime and daytime flux partitioning methods can yield similar or dissimilar daily RE, and sometimes even the seasonality can diverge significantly between them (e.g., US-PFa 2014 in Fig. 6).

Figure 8 Seasonal dynamics of GPP from daytime and nighttime partitioning method for 11 sites representative of different biomes.

Figure 9 The scatter plots of different phenology metrics from GPP data using daytime and nighttime partitioning methods.

The differences in flux seasonality metrics based on the partitioning method (Fig. 7) are sizable, non-systematic, and generally greater for season length metrics than transition date metrics (because the season length is determined by two transition dates, each of which is subject to deviation among methods). Interestingly, the seasonality metrics of flux development period (Fig. 7a, d, e, f) are more consistent than those of the peak flux period (Fig. 7g, h, i) and flux recession period (Fig. 7b, j, k, l). A more detailed analysis of the performance of the consistency of the different partitioning methods is the subject of future studies, and we may find answers more from information pertaining more to the accuracy of the flux estimates than the seasonality estimates.

Gross primary production (GPP)

The seasonality of GPP time series differed less between the two flux partitioning methods than did RE (Fig. 8). It is also obvious that the seasonalities of RE and GPP for the same sites differed significantly in terms of the seasonal timing, symmetry, peak duration, and other aspects. A more detailed assessment of these differences is the subject of a forthcoming analysis (Yang and Noormets, 2021).

As suggested by the extent of overlap in the flux estimates on Fig. 8, the seasonality metrics of GPP from the two partitioning methods were also more consistent compared to RE (Fig. 9). The scatter was distributed around the 1:1 line and with little bias. However, similar to RE, the metrics of the flux development period of GPP were also more consistent among the partitioning methods than those of peak flux period and flux recession period. Yet, the rate of flux development was more variable between the methods than the rate of flux recession (Fig. 9p, r, s, u). And like for ER, the phase duration metrics of GPP were more variable than the transition dates.

5 Data availability

This database can be obtained at the US Department of Energy's (DOE) Environmental Systems Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE, https://doi.org/10.15485/1602532; Yang and Noormets, 2020).

6 Significance

The Flux Seasonality Metrics Database (FSMD) presented here summarizes important latent features contained in the land surface biogeochemical fluxes and is likely to allow novel insights into the functioning of the biosphere and assist in the development and validation of novel functionality in Earth system models. Improving the predictive capabilities of ecosystem biogeochemistry models on interannual and decadal scales remains a challenge, and variability in the seasonality of different fluxes has been recognized as a key uncertainty. Importantly, the seasonality of GPP in models is often forced to match observations with arbitrary coefficients (Straube et al., 2018), as the divergence of LE and GPP seasonality is not captured in common LAI-driven models (Wu et al., 2017; Restrepo-Coupe et al., 2017). By developing process-specific seasonality references, explicit validation of these fluxes becomes possible. In addition, the ability to discern shifts in seasonality from those in flux capacity or vegetation structure can also be important in correcting the attributions of observed changes.

A standardized flux seasonality metrics dataset can also support other seasonality assessment tools, like near-surface (e.g., PhenoCam network; Richardson et al., 2012) and remote optical sensors (Broich et al., 2015; White and Nemani, 2006; Ganguly et al., 2010; Gamon et al., 2016). However, as these methods purport to infer GPP from the greenness information, they are vulnerable to the same lags between leaf development, LE, and GPP that have undermined the model development mentioned above. In addition, pixel heterogeneity and clouds can reduce the potential of remote sensing approaches.

The FSMD was designed to capture and depict the seasonality of different ecosystem processes. FSMD will be updated within 6 months of each major release of the FLUXNET database, making it possible to quantify the differences and similarities between different ecosystem processes in their responses to changes in climatic conditions. Some potential applications of this dataset have been mentioned before, and there are likely many others. We anticipate that the FSMD will stimulate new research in global change and Earth science disciplines where land–atmosphere exchange dynamics play a central role.