the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Reducing errors on estimates of the carbon uptake period based on time series of atmospheric CO_{2}

### Theertha Kariyathan

### Ana Bastos

### Julia Marshall

### Wouter Peters

### Pieter Tans

### Markus Reichstein

High-quality, long-time-series measurements of atmospheric greenhouse gases show interannual variability in the measured seasonal cycles. These changes can be analyzed to better understand the carbon cycle and the impact of climate drivers. However, nearly all discrete measurement records contain gaps and have noise due to the influence of local fluxes or synoptic variability. To facilitate analysis, filtering and curve-fitting techniques are often applied to these time series. Previous studies have recognized that there is an inherent uncertainty associated with this curve fitting, and the choice of a given mathematical method might introduce biases. Since uncertainties are seldom propagated to the metrics under study, this can lead to misinterpretation of the signal. In this study, we use an ensemble-based approach to quantify the uncertainty of the derived seasonal cycle metrics. We apply it to CO_{2} dry-air mole fraction time series from flask measurements in the Northern Hemisphere. We use this ensemble-based approach to analyze the carbon uptake period (CUP: the time of the year when the CO_{2} uptake is greater than the CO_{2} release): its onset, termination and duration. Previous studies have diagnosed CUP based on the dates on which the detrended, zero-centered seasonal cycle curve switches from positive to negative (the downward zero-crossing date, DZCD) and vice versa (upward zero-crossing date, UZCD). However, the UZCD is sensitive to the skewness of the CO_{2} seasonal cycle during the net carbon release period. Hence, we develop an alternative method proposed by Barlow et al. (2015) to estimate the onset and termination of the CUP based on a threshold defined in terms of the first derivative of the CO_{2} seasonal cycle. Using the ensemble approach, we arrive at a tighter constraint to the threshold by considering the annual uncertainty; we call this the ensemble of first derivative (EFD) method. Further, using the EFD approach and an additional curve-fitting algorithm, we show that (a) the uncertainty of the studied metrics is smaller using the EFD method than when approximated using the timing of the zero-crossing date (ZCD), and (b) the onset and termination dates derived with the EFD method provide more robust results, irrespective of the curve-fitting method applied to the data.

_{2}, Atmos. Meas. Tech., 16, 3299–3312, https://doi.org/10.5194/amt-16-3299-2023, 2023.

Ongoing in situ measurements of the atmospheric CO_{2} mixing ratio have revealed an increase in CO_{2} mole fraction in the atmosphere. The increase in atmospheric CO_{2} due to release of carbon from fossil fuel burning and land-use change is buffered by net CO_{2} uptake by the ocean and land biosphere (Keeling, 1960). Since then, many studies have used high-precision measurements of greenhouse gases at MLO (Mauna Loa) and other sites across the globe to better understand the role of CO_{2} in the global climate (e.g., Langenfelds et al., 2002; Keeling et al., 2017; Barlow et al., 2016). The analysis of such atmospheric time series helps to identify and isolate the long-term trends, inter-annual variability and seasonality of climatically important greenhouse gases (Thoning et al., 1989). However, these measurement records contain gaps and are influenced by local fluxes or synoptic-scale variability, which induces noise on the underlying climate signals. Hence the use of filtering and curve-fitting techniques to obtain smooth and continuous data has been an inevitable part of such studies (Trivett et al., 1989). The choice of mathematical method for data processing can, however, introduce biases that can result in misinterpretation of the signal (Nakazawa et al., 1997; Tans et al., 1989; Pickers and Manning, 2015; Barlow et al., 2015).

Curve-fitting methods are often used to preprocess atmospheric time series for analysis. Three examples are found in the commonly used software packages, HPspline (Bacastow et al., 1985), CCGCRV (Thoning et al., 1989) and STL (Cleveland et al., 1990). Each of these methods produces a gap-filled time series that contains the important features of the atmospheric record; however the resultant fitted curves vary significantly from each other owing to differences in their response to gaps and outliers in the original data. Pickers and Manning (2015) addressed the sensitivity of scientific conclusions to the curve-fitting method used, by repeating a scientific study (Piao et al., 2008) using two additional curve-fitting method. Both studies looked at changes in the CO_{2} seasonal cycle zero-crossing date (ZCD) for 10 midlatitude to high-latitude Northern Hemisphere stations. The re-analysis by Pickers and Manning (2015) found that the major conclusion of Piao et al. (2008) was robust but that inferences at individual stations depended on the curve-fitting method. This was corroborated by Barlow et al. (2015), who used a wavelet-based curve-fitting method to illustrate the sensitivity of various key aspects of the seasonal cycle of CO_{2} time series to the curve-fitting approach. Thus, the impact of bias introduced by data processing methods can vary based on the data set used and the type of analyses performed. Each method has its strengths and weaknesses; hence Pickers and Manning (2015) argued that data must be analyzed with multiple approaches to ensure that results are robust and free from bias. Despite this recommendation, studies that focus on metrics of time series such as the ZCD or seasonal cycle amplitude usually use a single curve-fitting method for analysis (e.g., Park et al., 2019; Piao et al., 2018), which can lead to differences in the conclusions that are drawn. An example is the disagreement in the direction of the trend of the CO_{2} seasonal cycle amplitude (SCA) at Alert, Canada between Chan and Wong (1990) and Keeling et al. (1996), as shown by Pickers and Manning (2015).

Metrics derived from CO_{2} time series such as the seasonal cycle peaks can be highly sensitive to data gaps and noise. This is especially true for metrics associated with the growing season onset at higher-latitude sites, where CO_{2} shows flat or multiple peaks in winter (Barlow et al., 2015). Hence, deriving other metrics like the timing of the carbon uptake period (CUP) from the seasonal cycle maximum results in less robust estimates. The CUP is defined as the time of the year during which the CO_{2} uptake is greater than the CO_{2} release. The onset and termination of the CUP are marked by the spring maximum and late-summer minimum of the seasonal cycle, respectively. However, the seasonal cycle at many observational sites is characterized by a flat peak or multiple peaks in winter, making it difficult to estimate the start of the CUP. To avoid this problem, previous studies have used the comparatively more unambiguous ZCD to approximate the timing and duration of the CUP (e.g., Piao et al., 2008). The ZCDs are the two dates in the seasonal cycle when the detrended CO_{2} curve crosses the zero line (an imaginary line passing through 0 ppm in the detrended CO_{2} seasonal cycle). Note that this period starts later than the seasonal spring maximum and ends later than the summer minimum; i.e., it is shifted compared to the CUP definition above. This approximation is thus based on the assumption that if the shape of the seasonal cycle does not change significantly, a change in the phase at one point (e.g., maximum) of the seasonal cycle can be traced as a relative phase change at other points (Barichivich et al., 2012). However, the shape of the seasonal cycle changes from year to year, and the CUP approximated using the ZCD may be erroneous (Barlow et al., 2015). This is illustrated in Fig. 1.

Barlow et al. (2015) show that using the time derivative of a time series can provide a more robust estimate of the key dates that define the CUP, compared to the conventional use of ZCDs. A threshold of this time derivative as a fraction of peak uptake in mid-summer was shown to be a robust metric to define both start (threshold 25 %) and end (threshold 0 %) of the CUP in their study. They used a synthetic data experiment applying a linear trend with substantial interannual variations in amplitude (±25 %) and CUP (±10 d) to a $\mathrm{d}{\mathrm{CO}}_{\mathrm{2}}/\mathrm{d}t$ time series, to show that in the absence of transport, their method can capture the prescribed linear trend of the CUP. We expand on that work here by additionally creating an ensemble of fitted time series using residual bootstrapping on a loess-fit. For each ensemble member we calculate the first derivative, allowing us to determine the timing of the various start, end, and peak moments in the CUP; its duration; and the individual uncertainty on each metric for each individual year in the time series. We call this the ensemble of first derivative method (EFD method). The EFD method accounts for the random and non-linear changes from year to year in the CO_{2} time series, allowing for a better handling of outlier years (in mean or uncertainty), which potentially improves trend analyses of seasonal cycle changes. We apply the EFD method to long time series and a set of stations covering the low latitudes, midlatitudes and high latitudes.

We first use the EFD method to confirm that the CO_{2} ZCDs are not the best proxy for determining the timing and duration of the CUP, also when the newly derived uncertainty is considered. We then demonstrate that the EFD method is independent of the skewness of the seasonal cycle, and we optimize the threshold for the CUP onset and termination based on the first derivative. The derived uncertainty also reveals that the robustness of various metrics is site-dependent, with high latitudes being sensitive to the seasonal cycle maximum (also found in Barlow et al., 2015) and low-latitude sites sensitive to the upward zero-crossing date (UZCD) of the CO_{2} seasonal cycle. We also tested whether the EFD method is sensitive to the specific curve-fitting method applied by fitting the data with the commonly used CCGCRV method, which is a frequency-domain-based filter, similar to the wavelet transform approach of Barlow et al. (2015). The measurements used in this study are presented in Sect. 2, and the EFD method is presented in Sect. 3. The results and the discussion on the findings can be found in Sects. 4 and 5 respectively, and Sect. 6 summarizes the findings of this study.

We use the discrete CO_{2} dry-air mole fraction from flask measurements from 10 observational sites of the NOAA/ESRL network (Dlugokencky et al., 2019, 2020), ranging from 19 to 82^{∘} N latitude. Table 1 lists the station names, station codes, their locations and the studied time period for each station (longer time records are available for MLO and NWR (Niwot Ridge), but these years have large data gaps of a year or more and hence are not considered for analysis). At these observational sites, air is sampled in glass flasks under background conditions; hence the dry-air mole fractions from the air samples are representative of the zonal mean atmospheric composition (Langenfelds et al., 2002). These air samples are collected weekly in pairs for quality control, and pairs with a difference of less than 0.5 ppm between the two samples are flagged as good-quality data (“good pairs”) (Dlugokencky et al., 2019, 2020). For our analysis, we use the mean value of each pair considered as “good pairs” and exclude low-quality measurements, which introduces irregular gaps in the data. The mean seasonal cycle of the higher-latitude stations (above 45^{∘} N latitude, i.e., SHM (Shemya Island), BRW (Barrow Atmospheric Baseline Observatory), ZEP (Ny-Ålesund, Svalbard, Norway and Sweden) and ALT (Alert)), is characterized by a broader maxima or multiple peaks in winter. Some lower-latitude stations like MLO, MID (Sand Island, Midway) and NWR have distinct seasonal cycles with clearly defined maxima, while others, like ASK (Assekrem, Algeria), AZR (Terceira Island, Azores) and WIS (Weizmann Institute of Science), have broader peaks (Fig. 2).

## 3.1 Loess fitting

The time series of CO_{2} can be described as the superposition of different modes of variability, acting at different frequencies. A standard approach to extract these modes of variability from the observations (*X*_{obs}(*t*)) is to define

where *X*_{trend}(*t*) is the low-frequency component of the data, which captures variability on multi-annual timescales; *X*_{seas}(*t*) represents the seasonal cycle, which can be expressed in terms of a series of harmonics; and *R*(*t*) captures the remaining variability (Cleveland et al., 1990).
The data used in this study are provided at approximately weekly time steps and include gaps, as described above. We fill gaps and estimate daily values by fitting a series of curves described in Eq. (1) and use residual bootstrapping (Kreiss and Lahiri, 2012) to generate an ensemble of 500 fitted curves consistent with the observational data for uncertainty estimation. Figure 3 describes the steps involved in the curve fitting and uncertainty estimation. Each step is described in detail below.

First, we separate the long-term trend and mean seasonal cycles (*X*_{trend}(*t*)+*X*_{seas}(*t*)) with a second-degree polynomial and four harmonic sinusoidal functions respectively (Bacastow et al., 1985). The remaining variability, *R*(*t*), is referred to as the residuals, which we verified to not show autocorrelation. We then fit a smooth curve to the residuals using the “loess” (local regression) method, which smooths the data, taking into account the gap lengths in the data. The “Caret” package (Kuhn, 2020) in R provides a method for optimizing the smoothing parameter for the “loess” regression using a mathematical method called *k*-fold cross validation. The optimization is based on five repetitions of 10-fold (*k*=10) cross-validation, where the subsamples are randomly sampled with restitution. The optimized smoothing parameters are then used to fit a smooth curve to the residuals (*R*(*t*)). The resulting smoothed residuals (*R*_{s}(*t*)), which contain the remaining variability, are added back to the other components (*X*_{trend}(*t*)+*X*_{seas}(*t*)). This produces a continuous and smooth data set that preserves short-term variations.

## 3.2 CCGCRV fitting

CCGCRV is a curve-fitting method developed by Kirk Thoning and Pieter Tans (Global Monitoring Laboratory (GML), NOAA) in the late 1980s. The method fits a combination of polynomials and annual harmonics to the data to approximate the long-term variation and seasonal cycle. The short-term and interannual variability are retained by filtering the residuals from the fit using a low-pass filter. A detailed description of the routines used for fitting the data and filtering of residual can be found in Thoning et al. (1989). In this study we use the C language version of CCGCRV, freely available at ftp://ftp.cmdl.noaa.gov/pub/john/ccgcrv/ (last access: 22 June 2022), for the curve fitting and finally obtaining a detrended time series. The values chosen for the input parameters were taken from Table 2 of Pickers and Manning (2015), who optimized them by fitting artificial data (short-term cut-off period *f*_{s}: 250 d; long-term cut-off period *f*_{l}: 1500 d; number of harmonic terms: 4; degree of polynomial function: 3).

## 3.3 Ensemble generation

Further, for uncertainty estimation, we generate 500 bootstrap samples from the curve-fitted data. For this, we calculate the difference between the smoothed data and the observational data, which gives the new set of residuals for generating bootstrap samples. These residuals are resampled (with replacement) and added to the initial fitted curve, producing a resampled time series. The resampled time series is processed as described in the preceding sections to obtain a continuous and smooth data set with daily values. The residual resampling and further processing are iterated 500 times to create an ensemble of 500 slightly different de-trended time series (bootstrap samples) which are all consistent with the observations (Fig. 3 shows these steps for loess fitting). The classical bootstrapping method (where the observations are resampled) cannot be applied directly to a time series data as the resampling step fails to replicate the time-dependent structure. Hence, we use residual bootstrapping where bootstrapping is applied to the residuals obtained from fitting a model to the raw data. The resampled time series thus show the same time dependence as the observational data but are produced from the fitted curve and a random component from the residual resampling.

The ensemble of fitted curves is used to constrain the uncertainty in seasonal cycle metrics estimates. If the estimated metrics differ largely across the bootstrap samples, it indicates that the metric estimate is influenced by the inherent uncertainty in extracting a definitive seasonal cycle, by curve-fitting the discrete data. Hence, interpreting these metrics without accounting for this uncertainty can be misleading.

## 3.4 Ensemble of first derivative (EFD) method

At high-latitude measuring stations, the CUP extends from the seasonal cycle maximum in spring to the seasonal cycle minimum in late summer (Barichivich et al., 2012), driven by CO_{2} uptake by ecosystems in the Northern Hemisphere. There is large uncertainty in associating the seasonal cycle maximum with the onset of the CUP, and the definition of the maximum is very sensitive to the curve-fitting method (Barichivich et al., 2012). The uncertainty in associating the timing of a maximum to the start of the CUP is larger than associating it with the ZCD, especially if the seasonal cycle is characterized by a fairly flat peak, or multiple peaks during the winter (Piao et al., 2008). Hence, previous studies have used the ZCD and their difference as proxies for the onset, termination and duration of the CUP, respectively. However, the period between the ZCD includes the CO_{2} release period that does not directly affect the CUP (Fig. 1). Therefore, we use the alternate method proposed by Barlow et al. (2015) to determine the timing and duration of the CUP from the first derivative of the mole fraction data, which more closely correspond to the spring maximum and the late-summer minimum times. We then estimate the uncertainty in the different CUP estimates using the spread of the ensemble members.

For each ensemble member, we calculate the first derivative of the time series as a proxy for the rate of CO_{2} uptake or release. The first derivative is at its minimum when CO_{2} uptake is most intense and reaches zero at the peak or trough of the seasonal cycle, i.e., when the sign of the integrated large-scale CO_{2} flux changes. However, a peak or a trough (as indicated by zero first derivative in Fig. 4) might not correspond to the spring maximum or late-summer minimum if the peak is flat or there are multiple peaks in winter. The timing of the CUP should be such that it closely corresponds to the timing of the spring maximum and late-summer minimum.

To determine the onset and termination of the CUP from CO_{2} mole fractions, we define a threshold, based on an ensemble of the first derivative of the time series. We define the threshold in terms of the first derivative of the CO_{2} dry-air mole fraction measurements analogous to Barlow et al. (2015). The first derivative can be seen as a proxy for the flux (not an exact correspondence, as the seasonal cycle at each site is affected by the atmospheric transport). The threshold is defined as *X* % of the first derivative minimum, and *X* is determined separately for the onset and termination of the CUP. The onset/termination of CUP is defined as the closest point to the threshold value before/after the first derivative minimum (Fig. 4). The threshold for the onset and termination is chosen such that (1) the uncertainty in the timing of onset and termination is minimized across the ensemble members, and (2) it represents as long a period as possible within the CUP. We varied the value of the parameter *X* until we find the optimum threshold. When *X* is 0 %, it corresponds to the time period between the seasonal cycle maximum and minimum, including the full CUP, but additional non-CUP periods may be erroneously included due to multiple peaks or flat maxima. By increasing the value of *X* we remove this error but can also truncate part of the “actual” CUP. Hence, we try to select a low value of *X* while reducing the uncertainty in the timing of the CUP.

For the EFD method, we first optimize the threshold as described in Sect. 3.4. By continuously increasing *X* we found the optimum for the termination is 0 %, and for the onset it is 12 %–13 %, with maximum CUP representation and no further reduction in the uncertainty beyond it. To be on the safe side, we chose 15 % as the threshold (for onset) in all our analyses. Incidentally, previous studies using flux measurements have also used 15 % of the maximum GPP as a threshold to define the start of the growing season (e.g., Wang et al., 2019). The result from varying *X* in steps between 0 %–20 % is shown in Fig. 5. When *X* is set to 0 %, the onset corresponds to the maximum of the seasonal cycle, hence the large spread in CUP onset at BRW and ALT. The interquartile range of 15.0 and 11.2 d respectively can be attributed to the multiple peaks or broad peak of the CO_{2} seasonal cycle at these stations.
Compared to the onset (30.8 d shown as whiskers of the box plots in Fig. 5a), the variability in the termination of the CUP is smaller, with a maximum range of 3.9 d (whiskers of the box plots in Fig. 5b). The standard deviation in the termination is highest at WIS where the median of the box plot for different threshold values is within 2.8 ± 0.2 d. Hence, to define the termination, we chose a threshold such that the standard deviation is minimized at WIS. This is achieved when *X* is set to 0 % (the median of the spread is then 2.5 d).

We estimate the duration of the CUP for each year using different approaches: (1) the difference between the seasonal cycle maximum and minimum times, (2) the difference in the ZCD and (3) using the EFD method. Figure 6a shows the duration of CUP for all the studied sites across the years (median of the 500 ensemble members), estimated using the three different methods. The size (interquartile range) of the box plots varies strongly across the stations for CUP duration calculated using the ZCD and the maxima and minima. At the lower-latitude stations MLO, MID and NWR, the variability in the CUP duration is larger than at the other stations when using the ZCD (“loess ZCD” in Fig. 6a). This is seen in the interquartile range of the “loess ZCD” box plots, with values of 43, 34 and 37 d for MLO, MID and NWR respectively, which is larger than for the other stations (within 15.2 ± 5.1 d). The large interquartile range of the CUP duration estimates using maximum and minimum times at the high-latitude sites BRW (46.75 d) and ALT (32 d) (“loess max.min” in Fig. 6a) mainly follows from the large variability in the timing of the seasonal cycle maximum across the ensemble members (Fig. 7).

When using the EFD method, the CUP estimates have least uncertainty across the ensemble members (Table 2). Figure 6b compares the standard deviation of the CUP duration across years at all studied sites and methods. The standard deviation is smaller when the EFD method is used for calculating the CUP duration, implying that this metric is better determined. The interquartile range of standard deviation is largest for the method using the dates of the seasonal cycle maximum and minimum, especially at higher-latitude stations like BRW (17 d) and ALT (13.7 d) and a lower-latitude station like WIS (17 d). At MLO, MID and NWR, using the ZCD to approximate the CUP duration results in a larger standard deviation (median of the spread is 5.5, 6.7 and 8.1 d respectively) in the CUP duration relative to the other methods used (the median of the spread for the other methods is within 1.75 ± 0.6 d).

Here we show that using the EFD method, the uncertainty in the CUP estimate is reduced across all the studied sites. Previous studies (Barichivich et al., 2012; Barlow et al., 2015) also noted the large uncertainty in using the seasonal cycle maximum and minimum to determine the CUP, which is similar to our result (Fig. 6b). Therefore this method will not be considered in further analysis here. However, when the ZCD is used to approximate CUP duration, there is also large uncertainty at the lower-latitude stations (e.g., the interquartile range for the “loess ZCD” box plot for MLO is 43 d, Fig. 6a). Nevertheless, the ZCD is a widely used approach (e.g., Piao et al., 2008; Barichivich et al., 2012, 2013), therefore the EFD method is compared against the ZCD here. The difference between the CUP estimates, using the two different methods (EFD and ZCD) varies from year to year, suggesting that the estimates cannot be corrected by simply adding an offset (Fig. 8). The *x*-axis range, showing the CUP from the ZCD in Fig. 8, has large year-to-year variation in the CUP, with the largest variation at MLO, NWR and ASK.

To further test the robustness of the CUP estimates based on the loess-fitted residual bootstrap method, we compared them against the CUP estimates from an ensemble using the CCGCRV curve-fitting method. Comparable results were obtained when the same CUP estimation method (ZCD/EFD) was applied to the ensemble members using the two different curve-fitting methods (Fig. 9a). The CUP duration calculated from the CCGCRV ensemble using the ZCD and the EFD method was within the range of their corresponding estimates from the loess-fitted ensemble members. The mean difference between the median of the “CCG ZCD” and “loess ZCD” estimates in Fig. 9a is 1.1 d, and between the median of “CCG EFD” and “loess EFD' estimates, the mean difference is only 0.6. However, the range of box plots corresponding to “loess ZCD” is larger relative to the “CCG ZCD”, resulting from the curve-fitting details. In the loess method, the long-term trend in the data is separated by fitting a quadratic polynomial, the decadal variability in the data is then retained which influences the ZCD, leading to more variability in the CUP approximated using the ZCD. The EFD method is less sensitive to the choice of the curve-fitting method (Fig. 9a) shown in the comparable “CCG EFD” and “loess EFD” numbers. Furthermore, we show that for both curve-fitting methods the standard deviation in the CUP duration estimate across the ensemble members is lowest for EFD method (Fig. 9b). Thus, the EFD method produces robust results irrespective of the particular curve-fitting method.

The CUP duration approximated using the ZCD shows larger spread for sites like MLO (with an interquartile range of 16 d for CCGCRV fitted data and 43 d for loess fitted data), irrespective of the curve-fitting method used. This is attributed to variability in the UZCD due to the skewness of the seasonal cycle during periods of net release and is similar in both the curve-fitting methods. Furthermore, we find that using the EFD method of CUP estimation resulted in smaller spread across the bootstrap samples for both the curve-fitting methods (Fig. 9b). This suggests that the period within the onset and termination defined by the EFD method, which includes only part of the drawdown period, is less variable than the time period between the ZCD, which includes parts of both the drawdown and release periods.

## 5.1 Uncertainty estimation with EFD method

In this study, we quantify the uncertainty in the CO_{2} seasonal cycle curve fitting by creating multiple residual bootstrap samples. The spread of this ensemble provides a measure of the uncertainty in the estimation of seasonal cycle metrics. The ensemble members are consistent with the observational data; hence we consider the variability of the metric estimate across the ensemble as a measure of uncertainty. The ensemble approach allows us to quantify the year-to-year change in different seasonal cycle metrics, and we see that the sensitivity of these metrics to curve fitting differs across latitudes and from year to year. Here we show that CO_{2} seasonal cycle metric estimates can be strongly sensitive to the method used; hence any method must be thoroughly evaluated before it can be used to derive trends from the atmospheric data. In Barlow et al. (2015) the robustness of the first derivative is tested by evaluating its ability to capture a known trend from a synthetic time series. They found a larger threshold value for the onset (25 %, suggesting a shorter CUP in their approach) from a synthetic data trend analysis in which they applied a linear trend with Gaussian variations of the peak uptake date to a CO_{2} time series. However, we argue that the data-derived year-to-year uncertainty from our ensemble provides a more robust threshold estimate and we derived a tighter threshold than Barlow et al. (2015) (15 % for onset). Further, Barlow et al. (2015) showed that their method can retrieve the true linear trend to within 10 %–25 %. Our EFD approach provides uncertainty on the year-to-year variability in the seasonal cycle metrics based on the fitted data residuals, which could be used in a trend analysis to give differential weights to each year. Also, trend analysis on the individual ensemble members would allow uncertainty on the trend to be calculated. Our demonstration of the EFD method on the CUP could be extended to other metrics that are derived directly from the seasonal cycle in a similar way, for example the peak-to-trough amplitude, especially when curve-fitting discrete data, at sites with broader or multiple peaks. In a similar fashion, the ensemble-based approach could be used to evaluate a newly proposed method or select an optimal method for evaluating any other metric based on reduced uncertainty.
The EFD approach, based on residual bootstrapping, assumes that the residuals are independent and as such may not be appropriate for hourly or daily in situ data. There is considerable auto-correlation between consecutive days in daily measurements, limiting the number of independent events to five or six per month, which is comparable to the scale of our weekly measurements. The uncertainties in curve fitting related to gaps in the data are also lower in the case of more frequent measurements.

## 5.2 Latitudinal dependence of metrics

The shape of the CO_{2} seasonal cycle varies with latitude. At the higher-latitude stations the seasonal cycle has a broader peak or multiple peaks in winter, and the timing of the seasonal cycle maximum cannot be interpreted as the onset of the CUP. Further, our analysis shows that there is large uncertainty in the timing of the seasonal cycle maximum (Fig. 7) at higher-latitude stations, which is in agreement with previous studies (e.g., Barichivich et al., 2012; Piao et al., 2008; Barlow et al., 2015). For example, for BRW shown in Fig. 10a, the atmospheric mixing ratios have a nearly constant value from January to May followed by a decrease in CO_{2} until a minimum is reached in August, also illustrated in Barlow et al. (2015) (their Fig. 6). If the seasonal cycle were determined solely by the biospheric fluxes, then the onset and termination of the CUP would be marked by the timing of the seasonal cycle maximum and minimum respectively (Barichivich et al., 2012). However, it can be noted that the estimated timing of the seasonal cycle maximum varies greatly across the bootstrap samples in BRW (inset of Fig. 10a), where the seasonal cycle is characterized by a flat peak. An earlier peak is likely to be associated with a flat maximum or multiple peaks that may result from transport (Parazoo et al., 2008) or other fluxes, rather than indicating the onset of the uptake period (Barlow et al., 2015). The timing of the ZCD at BRW is consistent across the ensemble members, which suggests that the timing of the ZCD (upward and downward) is less ambiguous. Other higher-latitude sites like ALT, SHM and ZEP and lower-latitude sites like ASK, WIS and AZR exhibit similar seasonal cycles, characterized by flat or multiple peaks and less ambiguous ZCD. However, the ZCDs are closer to the timing of the maximum uptake and release of CO_{2} (Manning, 1993) than to the actual onset and termination of the CUP. For example, in Fig. 1, the onset of the CUP occurs in June; however the downward zero crossing occurs in early July, and thus the CUP approximated using the ZCD explicitly excludes the start of the drawdown period.

At lower-latitude stations like MLO, MID and NWR, we find that the ZCD can vary across ensemble members as shown in Fig. 7, and in such cases the ZCDs are clearly not the best proxy for estimating the duration of the CUP. This is especially the case for the time series at MLO (Fig. 10b), which shows relatively a large spread of 5 d (median of the ensemble spread, rounded to the nearest integer) in the timing of the UZCD across the ensemble members. The seasonal cycle at MLO has well-defined peaks and troughs; hence the timings of the seasonal cycle maximum and minimum show only a small spread of 1 d (median of the ensemble spread, rounded to the nearest integer) across the bootstrap samples (inset of Fig. 10b). In this case, the EFD method gives a more robust estimate of the CUP duration.

We find that in addition to having a larger annual uncertainty, the range of CUP values over the study period for the ZCD approach is much larger than that of the EFD approach for some sites (Fig. 8). For example, at MLO the zero-crossing-approximated CUP ranges from 100 to 250 d, corresponding to a period of 3–8 months. Changes in the length of the growing season in the Northern Hemisphere are not expected to be this large. As an example, Jeong et al. (2011) estimated the length of the growing season using satellite measurements of the normalized difference vegetation index (NDVI). When integrating over the temperate Northern Hemisphere, the length of the phenology-derived growing season was found to vary by less than 25 d from 1982–2008. The ZCD approach includes changes in both the latter part of the net uptake period and the early-release period, making it difficult to separate the contribution of the net uptake and net release periods to the changes in the CUP estimate. To understand this large spread in CUP, we compare 2 years with very different CUP values estimated by the ZCD at MLO, 1992 with 192 d and 1998 with 147 d (Fig. 11). We find that the difference in the CUP estimate is due to the change in the early-release period, whereas the uptake periods are essentially the same. When using the EFD method, by contrast, the 2 years show similar CUP, 134 and 126 d, respectively. By definition, the EFD is not affected by differences in the net release period and therefore provides more robust CUP duration estimates.

## 5.3 Sensitivity of different methods to changes in the CUP

To test the sensitivity of the different methods to changes in the CUP, we performed a synthetic data experiment. We compare the EFD and ZCD methods and additionally an approach where we determine the CUP using a 25 % threshold as in Barlow et al. (2015) but using a single loess curve fitting rather than an ensemble. We manipulated the posterior net ecosystem exchange (NEE) flux by imposing pre-defined CUP changes (ΔCUP_NEE) for land pixels in the Northern Hemisphere that have a clearly defined seasonal cycle, using a randomly selected baseline year from the Jena Carboscope CO_{2} Inversion (Rödenbeck et al., 2003) (version ID: sEXT_ocNEET_v2021). The manipulated fluxes were then transported forward using an atmospheric transport model (TM3; Heimann and Körner, 2003). We then sampled the atmospheric mixing ratios at time steps corresponding to flask data measurements at the station locations. Then, we evaluated the CUP changes deduced from the simulated CO_{2} mixing ratio (ΔCUP_MR) in response to the imposed ΔCUP_NEE with the three different CUP estimation methods. The experiment shows that the EFD method is more sensitive to the changes in CUP_NEE compared to the ZCD method, which is seen by the higher slope of ΔCUP_MR to ΔCUP_NEE (0.312 vs. 0.08, respectively). Further, the uncertainty in the sensitivity of the ZCD and EFD methods was estimated by calculating the 95 % confidence intervals of the slopes of the regression lines by fitting multiple (10 000) sets of randomly resampled CUP_MR values (from its distribution) corresponding to each CUP_NEE value (shown in Fig. 12 for MLO). The confidence interval is large for the ZCD method, from −0.195 to 0.355. This implies that ΔCUP_NEE cannot be reliably inferred from ΔCUP_MR using the ZCD.
Further, we observed that the approach in Barlow et al. (2015) is similarly sensitive to changes in CUP_NEE as the ensemble-based approach; however the uncertainty in its sensitivity cannot be calculated as only a single CUP_MR value is available for every CUP_NEE, clearly showing the advantage of the ensemble approach to constrain uncertainties in metrics.

As seen from the experiment, the ZCD is the least sensitive to changes in the flux. Hence, years with extreme CUP approximated by the ZCD as in Fig. 8 may not indicate actual changes in the CUP but rather other factors, like reduced net respiration in the early-release period, which postpone the UZCD. This shift could be due to reduced net respiration in the early-release period, prolonging the time to reach the UZCD. This change is determined by the interplay of the CO_{2} uptake and release processes, which are influenced by physical factors like temperature, soil moisture and solar radiation. For example, in dry conditions there is less respiration by plants and slower decomposition of organic matter in the soil, resulting in reduced CO_{2} release to the atmosphere (Yan et al., 2018). The rate of decomposition further depends on the snow cover and available detritus content in the soil following leaf senescence. Furthermore, in the early-release period, when the solar radiation is not limiting, plants may continue to photosynthesize depending on water availability and temperature, leading to reduced net CO_{2} release. Thus, in years with extreme CUP as approximated by the ZCD, the physical processes that affect the release period should be investigated. In comparison to the CUP definition, the approximation by ZCD is also sensitive to variations after the summer minimum, i.e., during the early-release period.

The EFD and the approach in Barlow et al. (2015), which consider the first derivative of the concentration time series as a proxy for the large-scale spatially integrated flux (Barlow et al., 2015), are relatively more sensitive to the flux changes. However, the inferred changes cannot be directly translated to the underlying flux fields. The concentration time series do not exactly replicate the changes imposed in CUP_NEE, as these changes do not occur simultaneously across the Northern Hemisphere and hence do not represent a single changed flux but a superposition of many. These signals dampen each other when aggregated in CUP_MR, as seen in Fig. 12.

The synthetic experiment presented here used simulated data from only 1 year, thus controlling for the influence of inter-annual variability in atmospheric transport. In reality, atmospheric transport plays an important role in explaining a significant portion of observed CO_{2} variations at surface stations (e.g., Krol et al., 2018; Fu et al., 2015), which might affect any interpretation of the CUP metrics. An extensive study was carried out by Lintner et al. (2006), confirming the importance of atmospheric transport to account for some of the inter-annual variations in CO_{2} observed at Mauna Loa. Murayama et al. (2007) showed how year-to-year changes in the atmospheric transport create significant inter-annual variations in the downward zero-crossing date of the CO_{2} seasonal cycle that cannot be neglected. The ZCD is influenced by transport variability in both the late uptake and early-release periods. Hence, changes in the early-release period could be erroneously interpreted as changes in the CUP when using the ZCD. Also, atmospheric transport can contribute to the inter-annual variability in CUP estimates when using the EFD method. Hence, we recommend that the contribution of atmospheric transport at the studied background sites should be evaluated before interpreting and relating the CUP metrics to sources/sinks.

Here, we discuss a method for estimating the timing, duration and uncertainty of the CUP and related metrics from a discrete time series of CO_{2} dry-air mole fraction data. The uncertainty in the metrics is quantified using an ensemble of fitted time series generated through residual bootstrap sampling, a novel addition to the method presented in Barlow et al. (2015). Previous studies have used the timing of the ZCD as a proxy for defining the CUP; however the timing of the UZCD is influenced by the shape of the seasonal cycle, leading to large variability in the estimated CUP duration across the ensemble members for a given year for some of the studied sites, particularly at lower latitudes. The spread in the CUP duration across the ensemble members for a given year (i.e., the annual uncertainty) is lower for all studied sites when calculated using the EFD method. The EFD method depends directly on the timing and rate of the maximum CO_{2} uptake; hence the method is not affected by the shape change of the seasonal cycle outside the time period during which the CO_{2} uptake is larger than the CO_{2} release. Further, a synthetic data experiment showed that the approaches based on the first derivative of the CO_{2} mixing ratio are more sensitive to changes in CUP evaluated at MLO than the ZCD method. With the EFD method, the onset and termination are tightly constrained by considering the year-to-year change in the seasonal cycle. To test the impact of the curve-fitting method used, we generated bootstrap samples using both loess-fitted residuals and CCGCRV. The CUP duration estimated using the EFD method results in smaller spread for both curve-fitting methods. Further, for both curve-fitting methods, the standard deviation in the estimates across the ensemble members is smaller when using the EFD method, suggesting that the EFD method gives robust estimates. Thus, the EFD method allows for a robust estimate of the CUP that better reflects the CO_{2} drawdown period. This approach could be extended to other metrics of seasonal cycle analysis or to other curve-fitting methods, as was shown with the comparison to the CCGCRV results.

The code is made freely available under a Creative Commons-BY license, along with the documentation in this paper (https://doi.org/10.17617/3.ZKX9JS, Kariyathan, 2023). The CO_{2} dry-air mole fraction data for ALT, BRW and MLO (Dlugokencky et al., 2019) and for the other stations used in this study (Dlugokencky et al., 2020) are available from https://doi.org/10.15138/wkgj-f215.

The coding and analysis were performed by TK with contributions of MR. The study was conceptualized by JM, AB and MR with contributions from WP. JM, AB, WP, MR and PT contributed with expert knowledge. The original manuscript was drafted by TK, which was reviewed and edited by WP, JM, AB, MR and PT.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank NOAA/ESRL for the measurement and maintenance of the CO_{2} dry-air mole fraction data. We acknowledge the International Max Planck Research School for Global Biogeochemical Cycles (IMPRS-gBGC) for supporting the study.

The article processing charges for this open-access publication were covered by the Max Planck Society.

This paper was edited by Abhishek Chatterjee and reviewed by two anonymous referees.

Bacastow, R. B., Keeling, C. D., and Whorf, T. P.: Seasonal amplitude increase in atmospheric CO_{2} concentration at Mauna Loa, Hawaii, 1959–1982, J. Geophys. Res.-Atmos., 90, 10529–10540,
https://doi.org/10.1029/JD090iD06p10529, 1985. a, b

Barichivich, J., Briffa, K. R., Osborn, T. J., Melvin, T. M., and Caesar, J.: Thermal growing season and timing of biospheric carbon uptake across the Northern Hemisphere, Global Biogeochem. Cy., 26, GB4015, https://doi.org/10.1029/2012GB004312, 2012. a, b, c, d, e, f, g

Barichivich, J., Briffa, K. R., Myneni, R. B., Osborn, T. J., Melvin, T. M., Ciais, P., Piao, S., and Tucker, C.: Large-scale variations in the vegetation growing season and annual cycle of atmospheric CO_{2} at high northern latitudes from 1950 to 2011, Glob. Change Biol., 19, 3167–83, https://doi.org/10.1111/gcb.12283, 2013. a

Barlow, J. M., Palmer, P. I., Bruhwiler, L. M., and Tans, P.: Analysis of CO_{2} mole fraction data: first evidence of large-scale changes in CO_{2} uptake at high northern latitudes, Atmos. Chem. Phys., 15, 13739–13758, https://doi.org/10.5194/acp-15-13739-2015, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w

Barlow, J. M., Palmer, P. I., and Bruhwiler, L. M.: Increasing boreal wetland emissions inferred from reductions in atmospheric CH_{4} seasonal cycle, Atmos. Chem. Phys. Discuss. [preprint], https://doi.org/10.5194/acp-2016-752, in review, 2016. a

Chan, Y. H. and Wong, C. S.: Long-term changes in amplitudes of atmospheric CO_{2} concentrations at Ocean Station P and Alert, Canada, Tellus B, 42, 330–341, https://doi.org/10.1034/j.1600-0889.1990.t01-4-00003.x, 1990. a

Cleveland, R. B., Cleveland, W. S., McRae, J. E., and Terpenning, I.: STL: A Seasonal-Trend Decomposition Procedure Based on Loess (with Discussion), J. Off. Stat., 6, 3–73, 1990. a, b

Dlugokencky, E., Mund, J. W., Crotwell, A. M., Crotwell, M. J., and Thoning, K. W.: Atmospheric Carbon Dioxide Dry Air Mole Fractions from the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network, 1968–2018, Version: 2019-07, NOAA [data set], https://doi.org/10.15138/wkgj-f215, 2019. a, b, c, d, e, f

Dlugokencky, E., Mund, J. W., Crotwell, A. M., Crotwell, M. J., and Thoning, K. W.: Atmospheric Carbon Dioxide Dry Air Mole Fractions from the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network, 1968–2019, Version: 2020-07, NOAA [data set], https://doi.org/10.15138/wkgj-f215, 2020. a, b, c, d, e, f, g, h, i, j

Fu, Q., Lin, P., Solomon, S., and Hartmann, D. L.: Observational evidence of strengthening of the Brewer-Dobson circulation since 1980, J. Geophys. Res.-Atmos., 120, 10214–10228, https://doi.org/10.1002/2015JD023657, 2015. a

Heimann, H. and Körner, S.: The global atmospheric tracer model TM3, Technical Reports – Max-Planck-Institut für Biogeochemie, vol. 5, 131 pp., https://doi.org/10.4126/98-004424387, 2003. a

Jeong, S.-J., Ho, C.-H., Gim, H.-J., and Brown, M. E.: Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008, Glob. Change Biol., 17, 2385–2399, https://doi.org/10.1111/j.1365-2486.2011.02397.x, 2011. a

Kariyathan, T.: Reducing errors on estimates of the carbon uptake period based on time series of atmospheric CO_{2}, Edmond, V1 [code], https://doi.org/10.17617/3.ZKX9JS, 2023 a

Keeling, C. D.: The Concentration and Isotopic Abundances of Carbon Dioxide in the Atmosphere, Tellus, 12, 200–203, https://doi.org/10.1111/j.2153-3490.1960.tb01300.x, 1960. a

Keeling, C. D., Chin, J. F. S., and Whorf, T. P.: Increased activity of northern vegetation inferred from atmospheric CO_{2} measurements, Nature, 382, 146–149, 1996. a

Keeling, R. F., Graven, H. D., Welp, L. R., Resplandy, L., Bi, J., Piper, S. C., Sun, Y., Bollenbacher, A., and Meijer, H. A. J.: Atmospheric evidence for a global secular increase in carbon isotopic discrimination of land photosynthesis, P. Natl. Acad. Sci. USA, 114, 10361–10366, https://doi.org/10.1073/pnas.1619240114, 2017. a

Kreiss, J.-P. and Lahiri, S. N.: 1 – Bootstrap Methods for Time Series, in: Time Series Analysis: Methods and Applications, vol. 30 of Handbook of Statistics, edited by: Subba Rao, T., Subba Rao, S., and Rao, C., Elsevier, 3–26, https://doi.org/10.1016/B978-0-444-53858-1.00001-6, 2012. a

Krol, M., de Bruine, M., Killaars, L., Ouwersloot, H., Pozzer, A., Yin, Y., Chevallier, F., Bousquet, P., Patra, P., Belikov, D., Maksyutov, S., Dhomse, S., Feng, W., and Chipperfield, M. P.: Age of air as a diagnostic for transport timescales in global models, Geosci. Model Dev., 11, 3109–3130, https://doi.org/10.5194/gmd-11-3109-2018, 2018. a

Kuhn, M.: caret: Classification and Regression Training, r package version 6.0-85, Astrophysics Source Code Library, https://CRAN.R-project.org/package=caret (last access: 25 March 2020), 2020. a

Langenfelds, R. L., Francey, R. J., Pak, B. C., Steele, L. P., Lloyd, J. Trudinger, C. M., and Allison, C. E.: Interannual growth rate variations of atmospheric CO_{2} and its *δ*^{13}C, H_{2}, CH_{4}, and CO between 1992 and 1999 linked to biomass burning, Global Biogeochem. Cy., 16, 21-1–21-22, https://doi.org/10.1029/2001GB001466, 2002. a, b

Lintner, B. R., Buermann, W., Koven, C. D., and Fung, I. Y.: Seasonal circulation and Mauna Loa CO_{2} variability, J. Geophys. Res.-Atmos., 111, D13104, https://doi.org/10.1029/2005JD006535, 2006. a

Manning, M. R.: Seasonal Cycles in Atmospheric CO_{2} Concentrations, in: The Global Carbon Cycle, edited by: Heimann, M., Springer Berlin Heidelberg, Berlin, Heidelberg, 65–94, https://doi.org/10.1007/978-3-642-84608-3_3, 1993. a

Murayama, S., Higuchi, K., and Taguchi, S.: Influence of atmospheric transport on the inter-annual variation of the CO_{2} seasonal cycle downward zero-crossing, Geophys. Res. Lett., 34, L04811, https://doi.org/10.1029/2006GL028389, 2007. a

Nakazawa, T., Ishizawa, M., Higuchi, K., and Trivett, N. B. A.: Two curve fitting methods applied to CO_{2} flask data, Environmetrics, 8, 197–218,
https://doi.org/10.1002/(SICI)1099-095X(199705)8:3<197::AID-ENV248>3.0.CO;2-C, 1997. a

Parazoo, N. C., Denning, A. S., Kawa, S. R., Corbin, K. D., Lokupitiya, R. S., and Baker, I. T.: Mechanisms for synoptic variations of atmospheric CO_{2} in North America, South America and Europe, Atmos. Chem. Phys., 8, 7239–7254, https://doi.org/10.5194/acp-8-7239-2008, 2008. a

Park, T., Chen, C., Macias-Fauria, M., Tømmervik, H., Choi, S., Winkler, A., Bhatt, U. S., Walker, D. A., Piao, S., Brovkin, V., Nemani, R. R., and Myneni, R. B.: Changes in timing of seasonal peak photosynthetic activity in northern ecosystems, Glob. Change Biol., 25, 2382–2395, https://doi.org/10.1111/gcb.14638, 2019. a

Piao, S., Ciais, P., Friedlingstein, P., Peylin, P., Reichstein, M., Luyssaert, S., Margolis, H., Fang, J., Barr, A., Chen, A., Grelle, A., Hollinger, D., Laurila, T., Lindroth, A., Richardson, A., and Vesala, T.: Net carbon dioxide losses of northern ecosystems in response to autumn warming, Nature, 451, 49–52, https://doi.org/10.1038/nature06444, 2008. a, b, c, d, e, f

Piao, S., Liu, Z., Wang, Y., Ciais, P., Yao, Y., Peng, S., Chevallier, F.,
Friedlingstein, P., Janssens, I. A., Peñuelas, J., Sitch, S., and Wang, T.:
On the causes of trends in the seasonal amplitude of atmospheric CO_{2}, Glob. Change Biol., 24, 608–616, https://doi.org/10.1111/gcb.13909, 2018. a

Pickers, P. A. and Manning, A. C.: Investigating bias in the application of curve fitting programs to atmospheric time series, Atmos. Meas. Tech., 8, 1469–1489, https://doi.org/10.5194/amt-8-1469-2015, 2015. a, b, c, d, e, f

Rödenbeck, C., Houweling, S., Gloor, M., and Heimann, M.: CO_{2} flux history 1982–2001 inferred from atmospheric data using a global inversion of atmospheric transport, Atmos. Chem. Phys., 3, 1919–1964, https://doi.org/10.5194/acp-3-1919-2003, 2003. a

Tans, P. P. K. W. T., Elliott, W., and Conway, T. J.: Background Atmospheric CO_{2} patterns from weekly flask samples at Barrow, Alaska: Optimal signal recovery and error estimates, in The Statistical Treatment of CO_{2} Data Records, NOAA Technical Memorandum, 173, 131, 112–123, https://www.arl.noaa.gov/documents/reports/arl-173.pdf (last access: 25 June 2023), 1989. a

Thoning, K. W., Tans, P. P., and Komhyr, W. D.: Atmospheric carbon dioxide at Mauna Loa Observatory: 2. Analysis of the NOAA GMCC data, 1974–1985, J. Geophys. Res.-Atmos., 94, 8549–8565, https://doi.org/10.1029/JD094iD06p08549, 1989. a, b, c

Trivett, N. B. A., Higuchi, K., and Symington, S.: Trends and seasonal cycles of atmospheric CO_{2} over Alert, Sable Island, and Cape St. James, as analyzed by forward stepwise regression technique, NOAA Technical Memorandum ERL ARL- 173, Air Resources Laboratory, Silver Spring, Maryland, USA, 173, 131, 27–42, https://www.arl.noaa.gov/documents/reports/arl-173.pdf (last access: 25 June 2023), 1989. a

Wang, X., Xiao, J., Li, X., Cheng, G., Ma, M., Zhu, G., Altaf Arain, M., Andrew Black, T., and Jassal, R. S.: No trends in spring and autumn phenology during the global warming hiatus, Nat. Commun., 10, 2389, https://doi.org/10.1038/s41467-019-10235-8, 2019. a

Yan, Z., Bond-Lamberty, B., Todd-Brown, K. E., Bailey, V. L., Li, S., Liu, C., and Liu, C.: A moisture function of soil heterotrophic respiration that incorporates microscale processes, Nat. Commun., 9, 2562, https://doi.org/10.1038/s41467-018-04971-6, 2018. a

_{2}measurements with noise and gaps. The CUP metrics derived with the approach are more robust and have less uncertainty than when estimated with the conventional methods.