Investigating the differences in calculating global mean surface CO₂ abundance: the impact of analysis methodologies and site selection

The World Meteorological Organization (WMO) Global Atmosphere Watch (GAW) coordinates high-quality atmospheric greenhouse gas observations globally and provides these observations through the WMO World Data Centre for Greenhouse Gases (WDCGG) supported by Japan Meteorological Agency. The WDCGG and the National Oceanic and Atmospheric Administration (NOAA) analyse these measurements using different methodologies and site selection to calculate global annual mean surface CO₂ and its growth rate as a headline climate indicator. This study introduces a third hybrid method named GFIT, which serves as an independent validation and open-source alternative to the methods described by NOAA and WDCGG. We apply GFIT to incorporate observations from most WMO GAW stations and 3D modelled CO₂ fields from CarbonTracker Europe (CTE). We find that different observational networks (i.e. NOAA, GAW, and CTE networks) and analysis methods result in differences in the calculated global surface CO₂ mole fractions equivalent to the current atmospheric growth rate over a 3-month period. However, the CO₂ growth rate derived from these networks and the CTE model output shows good agreement. Over the long-term period (40 years), both networks with and without continental sites exhibit the same trend in the growth rate (0.030 ± 0.002 ppm yr⁻¹ each year). However, a clear difference emerges in the short-term (1-month) change in the growth rate. The network that includes continental sites improves the early detection of changes in biogenic emissions.

1 Introduction

Global mean surface temperature averaged over 2011–2020 has increased by about 1.09 ^∘C relative to the average temperature of 1850–1900 (Gulev et al., 2021). The increasing amount of atmospheric carbon dioxide (CO₂), together with increases in other greenhouse gases, is the main driver of the warming (Eyring et al., 2021). After being relatively stable between 180 ppm (ice age) and 280 ppm (interglacial) for the past 800 000 years (Lüthi et al., 2008), the annual average CO₂ level of the atmosphere has increased since the industrial revolution from roughly 277 ppm in 1750 to 415.7 ± 0.2 ppm in 2021 (WMO, 2022), due to emissions of CO₂ related to human activities like burning of fossil fuels and land use changes (Friedlingstein et al., 2022). Mean global atmospheric CO₂ annual growth rate (G_ATM) is an important constraint on the global carbon cycle. Based on the most recent global carbon budget (GCB) analysis (Friedlingstein et al., 2022), the total emission of CO₂ due to human activities was 10.2 ± 0.8 GtC yr⁻¹ in 2020, of which 3.0 ± 0.4 GtC yr⁻¹ was captured by the ocean sink and 2.9 ± 1 GtC yr⁻¹ by the terrestrial sink, leaving a net increase of 5.0  ± 0.2 GtC yr⁻¹ of CO₂ in the atmosphere, corresponding to an atmospheric CO₂ mole fraction increase of 2.4 ± 0.1 ppm yr⁻¹. (The conversion factor comes from Ballantyne et al., 2012.)

As the atmosphere mixes the contributions of all sources and sinks, an observational global average CO₂ mole fraction can be constructed if there are enough observations to represent the spatial and temporal variations across the globe. Since most land masses are concentrated in the Northern Hemisphere, and the highest anthropogenic emissions (e.g. during winter) occur in the relatively narrow latitudinal band between 30 and 60^∘ N, relatively large spatial and temporal gradients in CO₂ mole fraction exist in and around that region. Due to convective and advective mixing, the average mixing time of air within the same latitudinal bands varies from several weeks to a month. However, mixing between latitudinal bands is slower, especially the exchange between the Northern and Southern hemispheres, which has an approximate interhemispheric transport time of 1.4 ± 0.2 years (Patra et al., 2011). The interplay of the latitudinal and interhemispheric differences in fossil fuel emissions and seasonal exchange with land biota (Denning et al., 1995) creates a latitudinal and interhemispheric gradient that requires a sufficiently dense network to capture a representative global annual mean.

However, measurement stations that are close to sources or sinks may not be representative of a large atmospheric volume and the average signal at their latitude. Therefore, inclusion of these observations might introduce biases on the global mean CO₂ and its growth rate. These biases can be avoided by filtering of data and a careful selection of spatially representative stations, as done by NOAA in their use of 43 stations (Fig. 1) that are considered to be representative for the Marine Boundary Layer (MBL reference network, https://www.esrl.noaa.gov/gmd/ccgg/mbl/mbl.html;last access: 7 December 2023). An additional data processing step developed by NOAA to further avoid biases due to unrepresentative local signals is filtering and smoothing, by using a combination of a low pass filter and decomposition into a fitted long-term trend and seasonal cycle (Thoning et al., 1989), hereafter referred to as the NOAA analysis. These fits can also be used to fill gaps for missing data, though care must be taken to avoid extrapolation errors before and beyond the time covered by the data record of the station. The WMO Global Atmosphere Watch (GAW) World Data Centre for Greenhouse Gases (WDCGG) publishes global averages mole fraction for CO₂ and the other major greenhouse gases in the annual WMO GAW Greenhouse Gas Bulletin (latest version: WMO, 2022). They use curve fitting and filter methods that are very similar to those developed by NOAA, but WDCGG includes continental locations that are potentially more influenced by local sources and sinks (Tsutsumi et al., 2009).

The NOAA MBL observations are all part of the NOAA cooperative global air sampling network and analysed in the same laboratory. All NOAA flask–air observations are traceable to the current scale WMO–CO₂–X2019 (Hall et al., 2021) that is maintained by NOAA Global Monitoring Laboratory (GML). In contrast, the WDCGG data originate from multiple independent laboratories (including NOAA GML), that together form a network of hundreds of stations coordinated by WMO GAW (http://gawsis.meteoswiss.ch; last access: 7 December 2023). Having a multitude of independent laboratories carries an additional risk of biases due to differences in sampling, measurement, and analysis methods, for example calibration scales, although much care is taken to avoid these by coordination in the network and use of a common calibration scale from WMO Central Calibration Laboratories (CCL) guided by a set of strict measurement compatibility goals (WMO, 2022). The different selection of stations results in a larger seasonal cycle amplitude in WDCGG results compared with those of NOAA and a small but quite consistent bias in global surface annual mean CO₂ mole fraction (Tsutsumi et al., 2009). The NOAA estimate of global surface annual mean CO₂ mole fraction is expected to be lower (e.g. ∼ 0.35 ppm lower than the WDCGG estimate, Tsutsumi et al., 2009) compared with a full global surface average because areas with large sources are not represented. However, the two aforementioned approaches neither represent the parts of the atmosphere with low CO₂ mole fraction levels (i.e. the full troposphere, up to ∼ 8–15 km altitude, and the stratosphere), nor do they cover the regions of the world with substantial observational gaps.

In this paper, we propose a data integration method to estimate the global mean surface CO₂ and its growth rate, named GFIT. This method serves as an independent validation of the methods as described by NOAA and WDCGG through a completely independent and open-source implementation. The global mean surface CO₂ refers to the mean CO₂ mole fraction within the planetary boundary layer, which extends from the earth's surface up to a few hundred or thousand metres in height. We apply the GFIT methodology to incorporate CO₂ data from the GAW network (139 stations; Fig. 1) and the modelled CO₂ distribution from a well-established 3D global transport model (TM5: Transport Model 5; Krol et al., 2005; Peters et al., 2004). We investigate the influence of small differences between the three methodologies and whether these are significant or not for calculating the global mean surface CO₂ and its growth rate, how consistent the GFIT and WDCGG approaches are with each other, and how they compare with NOAA analysis and estimates derived from a CO₂ simulation with the 3D transport model TM5. These 3D CO₂ results for 2001–2020 using TM5 are performed in the CarbonTracker Europe framework (CTE; Peters et al., 2004; Van Der Laan-Luijkx et al., 2017), where the CO₂ uptake and emission fluxes are optimized by the inversion system to minimize the mismatch between the in situ observations and the modelled CO₂ mole fraction. CTE generally has a good representation of the CO₂ field, with mean biases with respect to independent aircraft measurements of generally less than 0.5 ppm (Friedlingstein et al., 2022). Furthermore, the inferred CO₂ fluxes from CTE fit well within the ensemble of those of other inversions used for the evaluation of global carbon budget (e.g. Friedlingstein et al., 2022).

2 Methods and data

2.1 The WMO GAW observations and WDCGG analysis method

The WMO GAW network measurements are archived and distributed by WDCGG, hosted by the Japan Meteorological Agency. The GAW observations used in this study originate from 139 selected stations of the GAW network, and all observations are on the WMO standard scale WMO–CO₂–X2019. The details on the station selection are described in Tsutsumi et al. (2009), which mainly excludes stations located in the Northern Hemisphere that show large standard deviations from the latitudinal fitted curve. The remaining 139 stations show a more reasonable latitudinal scatter range (Fig. 1).

Figure 1 Three observation networks are employed to assess the impact of continental site inclusion when calculating global CO₂ mole fraction and its growth rates. The NOAA network (43 sites; yellow stars) comprises MBL sites only. The selected GAW global network (139 sites; red dots) includes both MBL sites and continental sites, for example from the Advanced Global Atmospheric Gases Experiment (AGAGE) and the European Integrated Carbon Observation System (ICOS) contribution network. The CTE network serves as the global network for the CTE model evaluations (230 sites; blue dots), and comprises MBL sites and a more extensive inclusion of continental sites.

The WDCGG global analysis method (hereafter WDCGG method), as described in Tsutsumi et al. (2009), includes the mentioned station selection, a data fitting and filter (involves data interpolation and extrapolation), and calculation of the zonal and global mean mole fractions, trends, and growth rates. The procedure is also summarized in Sect. S1 in the Supplement. The output from the global analysis by the WDCGG method is used to compare with an alternative method (GFIT) that we designed to follow as closely as possible the fit and filter method (Conway et al., 1994) deployed by NOAA and is described in the Sect. 2.3.

2.2 CTE model output and station observations

CarbonTracker Europe is a global model of atmospheric CO₂ and designed to keep track of CO₂ uptake and release at the earth's surface over time (Van Der Laan-Luijkx et al., 2017). CTE incorporates an off-line atmospheric transport module (TM5, Peters et al., 2004; Krol et al., 2005) driven by ECMWF ERA5 data, and there are four prescribed fluxes (i.e. from ocean, biosphere, fire, and fossil fuel), which are transported in the model, together with the transported initial CO₂ field. CTE also includes a data assimilation system that applies an ensemble Kalman filter to optimize the biogenic and ocean fluxes for a combination of plant-functional types and climate zones to improve the fit of the simulated concentrations with observations. The optimized fluxes from the data assimilation have been used in Global Carbon Project (GCP) 2021, and the comparison of CTE CO₂ product with the other data assimilation systems used in GCP shows good agreement (within 0.8 ppm at all latitude bands) (Friedlingstein et al., 2022).

The CTE model data used here consist of simulated monthly CO₂ mole fraction at 1 × 1^∘ horizontal resolution and 25 levels in the vertical direction, and the data period ranges from 2001 to 2020 which has no influence of model spin-up (Krol et al., 2018). From the CTE output a set of simulated synthetic atmospheric CO₂ mole fractions with monthly resolution can be extracted within grid cells where stations are situated. This study analyses monthly observation data (1980–2020) and synthetic time series (2001–2020) by using the GFIT method (Sect. 2.3) and attempts to estimate global mean CO₂ mole fraction and its growth rate. The observed CO₂ mole fractions are taken from 230 out of 290 global-wide distributed stations (Fig. 1; the station selection is summarized in Sect. S2). The data come from the GLOBALVIEW-plus V8 ObsPack data product (Schuldt et al., 2022) and include surface-based, shipboard-based, and tower-based measurements.

2.3 The GFIT method

The temporal pattern of CO₂ measurement records at locations around the globe can be explained as the combination of roughly three components: a long-term trend, a non-sinusoidal yearly cycle (or seasonality), and short-term variations. This study synchronizes monthly CO₂ records with the fitting and filter method developed at the NOAA Global Monitoring Laboratory (Conway et al., 1994; Thoning et al., 1989), without extrapolation. The station selection and CO₂ averaging method are kept the same as in the WDCGG method (Sect. S1). This method will be referred to as the GFIT method and will be compared with the WDCGG method without extrapolation. The only difference from the WDCGG method without extrapolation is the fitting and filter method. All code for the method described here was developed in Python and is available as a Jupyter notebook under a GPL license (https://doi.org/10.18160/Q788-9081, Wu, 2023). The GFIT method can be summarized and illustrated by the three steps described in the next subsections.

2.3.1 Fitting and filter

CO₂ records from each station can be abstracted as a combination of long-term trend and seasonality, which can be fitted by a function consisting of polynomial and harmonics. We applied a linear regression analysis based on three polynomial coefficients and four harmonics (Eq. 1) to fit CO₂ data using general linear least-squares fit (LFIT, Press et al., 1988):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}f\left(x\right)& =a_{\mathrm{0}}+a_{\mathrm{1}}t+a_{\mathrm{2}}{t}^{\mathrm{2}}+\mathrm{\dots }+a_k{t}^k\\ & +\sum _{n=\mathrm{1}}^{n_{\mathrm{h}}}\left(A_n\cos \mathrm{2}\mathit{\pi }nt+B_n\sin \mathrm{2}\mathit{\pi }nt\right),\end{array}\end{array}$$

where a_k, A_n, and B_n are fitted parameters; t is the time from the beginning of the observation and it is in months and expressed as a decimal of its year; k denotes the polynomial number, k=2; n_h denotes harmonic number; and n_h=4. Figure 2 illustrates the function fit to CO₂ data to obtain the annual oscillation (red line in Fig. 2a), is a combination of a polynomial fit to the trend (blue line in Fig. 2a), and is a harmonic fit to the seasonality (green line in Fig. 2b).

The residuals are the difference between raw data and the function fit (black dots in Fig. 2c). The filtering method is based on Thoning et al. (1989) which transforms CO₂ data from time domain to frequency domain using a fast Fourier transform (FFT), then applies a low pass filter to the frequency data to remove high-frequency variations, and then transforms the filtered data back to the time domain using an inverse FFT. The short-term (a cut-off value of 80 d; red line in Fig. 2c) and long-term (a cut-off value of 667 d; blue line in Fig. 2c) filters used here are the same as in NOAA method and applied to obtain the short-term and interannual variations that are not determined by the fit function. The original code is also available as Python code from the NOAA website (https://gml.noaa.gov/aftp/user/thoning/ccgcrv/; last access: 7 December 2023).

2.3.2 Calculate smoothed CO₂ and long-term trend

The result of filtering residuals is added to the fitted curve to obtain smoothed CO₂ and its long-term trend. The smoothed CO₂ comprises fitted trend, fitted seasonality and smoothed residuals (red line in Fig. 2d), the latter removes only short-term variations or noise. The long-term trend comprises fitted trend and residual trend, which removes seasonal cycle and noise (blue line in Fig. 2d).

2.3.3 Calculate CO₂ growth rate, G_ATM

The CO₂ growth rate (G_ATM) is determined by taking the first derivative of the long-term trend. However, the growth is made up of discrete points, e.g. the black dots in Fig. 3a show the trend points. In this case, a cubic spline interpolation is applied to the trend points, in which the spline curve passes through each trend point, such as the blue line in Fig. 3a. G_ATM is obtained by taking the derivative of the spline at each trend point (Fig. 3b).

Figure 2 Example of analysed CO₂ data from station Pallas (PAL, Finland), illustrating GFIT curve fitting and filter method. Panel (a) shows monthly averaged CO₂ (dots), curve fitting with 2-degree polynomial and 4-degree harmonics (red line), and long-term trend estimated by a 2-degree polynomial (blue line). Panel (b) shows seasonality estimated by 4-degree harmonics. Panel (c) shows the residuals of raw data from the function fit (black dots). The red line is obtained by the short-term filter and the blue line is obtained by the long-term filter. The cyan dots show the residuals of raw data from the sum of fitted curve and smoothed residuals. Panel (d) shows final processed CO₂, which comprises fitted trend, fitted seasonality, and smoothed residuals (red line). The blue line shows the final trend which comprises fitted trend and residuals trend.

Figure 3 Example of CO₂ growth rate. The raw data is the same as in Fig. 2 from station Pallas (PAL, Finland). Panel (a) shows the trend points (black dots) and its cubic spline interpolation (blue line). Only if there is missing data the black dots will be less, and the blue line is clearly visible. Panel (b) shows the G_ATM at each trend point.

3 Results

Global averaged surface CO₂ and its G_ATM are calculated using the WDCGG method and our GFIT method based on the data from the GAW and CTE networks (Fig. 1). The different observation networks and their analysis methods are listed in Table 1. We calculated the global means and its G_ATM by area-weighted averaging the zonal means over each latitudinal band (30^∘), following the same CO₂ averaging method as described in Tsutsumi et al. (2009). A bootstrap method is used to estimate the uncertainties of global CO₂ mean and its G_ATM, which is an almost identical uncertainty analysis as presented by Conway et al. (1994) who constructed 100 bootstrap networks for the NOAA analysis. We construct 200 bootstrap networks, consistent with the WDCGG analysis in Tsutsumi et al. (2009). For each bootstrap network, we randomly draw the same number of sites as the actual network (e.g. 139 sites for GAW network) with replacement from the actual network, which means some sites are missing whereas others are represented twice or more often. We calculate global mean CO₂ mole fraction and its G_ATM for each network and then calculate the statistics (i.e. mean and 68 % confidence interval (CI)) on the 200 networks. All uncertainties in this paper are reported as ± 68 % CI.

Table 1 Description of the three observation networks and their analysis methods.

3.1 Globally averaged surface CO₂ mole fraction and its G_ATM

Figure 4 presents a monthly comparison of globally and locally averaged CO₂ mole fractions and their G_ATM from 1980 to 2020. The statistical metrics assessing the agreement of these monthly comparisons are available in Fig. 5 (for 2001–2020) and Fig. S1 in the Supplement (for 1980–2020). The statistical metrics for the annual comparisons can be found in Fig. S2 (for 2001–2020) and Fig. S3 (for 1980–2020). They exhibit a similar pattern to the monthly comparisons (i.e. Figs. 5 and S1).

Figure 4 Comparison of globally and locally averaged CO₂ mole fraction (a) and its G_ATM (b) from 1980 to 2020. Panel (a) shows the global monthly CO₂ mole fraction from 139 GAW sites (estimated from observations only), 43 NOAA MBL sites, and those from 230 sites used in CTE (either from observations or model output). The two local CO₂ mole fractions are from Mauna Loa (MLO; cyan line) and South Pole (SPO; magenta line) stations, analysed using the GFIT method. The red and blue lines show the CO₂ derived from GAW (GFIT) and GAW (WDCGG), respectively. The green and orange lines show the CO₂ derived from CTE_obs (GFIT) and CTE_output (GFIT), respectively. The right y axis shows their difference from NOAA CO₂ mole fraction, and the dashed lines show the mean of the difference over the available period. Panel (b) compares the corresponding global and local CO₂ growth rate; the legend refers to (a). The shadow area shows the uncertainty as a 68 % confidence interval obtained by the bootstrap analysis.

Globally averaged monthly surface CO₂ mole fractions, derived from the GAW network (GAW (GFIT) or GAW (WDCGG)), are significantly (p<0.05) higher by 0.329–0.335 ppm during 1980–2020 (Fig. S1a) and 0.370–0.390 ppm during 2001–2020 (Fig. 5a) when compared with the NOAA analysis (Fig. 4a). This finding aligns with that of Tsutsumi et al. (2009), who reported a 0.350 ppm higher global average in the GAW network during 1983–2006. The higher estimate from the GAW network can be attributed to the inclusion of more diverse sites, encompassing not only NOAA's MBL sites but also additional continental sites (Fig. 1).

Both global CO₂ and its G_ATM derived from the GAW (GFIT) and GAW (WDCGG) are nearly overlapping (the red and blue lines) in Fig. 4a and b. The statistical metrics (Figs. 5 and S1) indicate a high agreement (ME < 0.020 ppm, RMSE < 0.145 ppm, and r>0.999 for CO₂ mole fraction; ME < 0.005 ppm yr⁻¹, RMSE < 0.108 ppm yr⁻¹, and r>0.982 for G_ATM) between these two methods, which confirms that the GFIT method agrees well with WDCGG method without extrapolation. The WDCGG method with extrapolation (i.e. GAW (WDCGG+)), which involves extrapolating the long-term trend of each station to match the period of the most long-running station and adding it to the average seasonal variation to synchronize data period of all stations (Tsutsumi et al., 2009), produces 0.096 ppm significantly (p<0.05) higher values than the global monthly surface CO₂ mole fraction derived from the GAW (WDCGG) during the common period 1984–2020 (Table S1 in the Supplement). However, the extrapolation has a minimal effect (RMSE = 0.076 ppm yr⁻¹ and ME = −0.011 ppm yr⁻¹; Table S1) on the CO₂ growth rate.

Globally averaged monthly surface CO₂ derived from CTE_obs (GFIT) and CTE_output (GFIT) are 0.422 ppm (1980–2020; Fig. S1) and 0.668 ppm (2001–2020; Fig. 5) significantly (p<0.05) higher compared with the NOAA analysis, respectively (Fig. 4a). Comparing the global mean of CTE_obs (GFIT) with CTE_output (GFIT) during the common period of 2001–2020, we observe a low bias (0.069 ppm in CTE_output; Fig. 5a), which suggests that the CTE model results can reasonably reproduce the global mean CO₂ levels. The global annual CO₂ mole fraction from CTE_obs (GFIT), CTE_output (GFIT), and CTE_global (GFIT) is 0.367, 0.299, and 0.186 ppm significantly (p<0.05) higher than the result of the GAW (GFIT), respectively (Fig. 5a). The higher global mean from CTE_obs (GFIT) and CTE_output (GFIT) can be attributed to the presence of more sites in the Northern Hemisphere within the CTE network compared with the GAW network. The lower bias observed between GAW (GFIT) and CTE_global (GFIT) suggests that the GAW network provides a good representation of the low-level atmosphere (i.e. 0 to 0.35 km altitude) at global scale, or the CTE model performs well in the low-level atmosphere.

A common approach to estimate global surface CO₂ mole fraction is by using one or two representative sites, such as MLO and SPO. The globally averaged monthly surface CO₂ mole fractions, derived from the GAW, CTE, and NOAA networks, are significantly (p<0.05) lower by 0.46–0.88 ppm during 1980–2020 (Fig. S1a) and 0.45–1.19 during 2001–2020 (Fig. 5a) than the local CO₂ estimates solely based on MLO measurements. Conversely, these global monthly CO₂ mole fractions are significantly (p<0.05) higher by 1.91–2.24 ppm during 1980–2020 (Fig. S1a) and 2.21–2.94 ppm during 2001–2020 (Fig. 5a) when compared with local measurements at the SPO site. Furthermore, the global seasonal cycle leads the local cycle at MLO by approximately 1 month (estimated by averaging the time difference between the peaks of their seasonal cycles). In contrast, the local cycle at SPO is not evident and is opposite to the global seasonal cycle (Fig. 4a).

Figure 5 Pair-wise statistical metrics assess the agreement of monthly global and local CO₂ mole fraction (ppm) and its G_ATM (ppm yr⁻¹) across various networks and methodologies (see Fig. 4 and Table 1) for the period 2001–2020. Panel (a) presents the mean error (ME) quantifying the difference for each pair, focusing on CO₂ mole fraction, while (b) does the same for G_ATM. The significance levels of paired t test for ME are indicated as follows: ^* p<0.1, ${}^{**}$ p<0.05 and ${}^{***}$ p<0.01. Panels (c) and (d) present the root mean squared error (RMSE) for CO₂ mole fraction and G_ATM, respectively. Panels (e) and (f) present the Pearson correlation coefficient (r) for CO₂ mole fraction and G_ATM, respectively.

Figure 6 Trend analysis of the global CO₂ growth rate from 1980 to 2020. Panel (a) shows the trends in CO₂ growth rate for the GAW network (red line), the CTE network (green line), and the NOAA network (black line) during the whole period 1980–2020. The CO₂ growth rate is derived from GAW (GFIT), CTE_obs (GFIT), and NOAA analysis (Fig. 4b). Panel (b) shows the trend in CO₂ growth rate for each month during 1980–2020, calculated as the derivative of the growth rate. The grey bands mark the period of three strong El Niño events, i.e 1987–1988, 1997–1998, and 2014–2016.

Despite differences in the global averaged surface CO₂ mole fractions derived from different networks and analysis methods, the G_ATM derived from GAW network, CTE network and its model output, and NOAA network exhibits strong agreement during 1980–2020 (ME < 0.031 ppm yr⁻¹, RMSE < 0.217 ppm yr⁻¹, and r>0.948; Figs. 4b and S1). The differences in the G_ATM remain below 0.023 ppm yr⁻¹ during 2001–2020, with low or no significance level (Fig. 5b), especially when comparing the annual G_ATM (Fig. S2b). Furthermore, over the long-term period of 40 years, the estimated local growth rate at MLO (ME < 0.046 ppm yr⁻¹ higher, RMSE < 0.272 ppm yr⁻¹, and r>0.915) and SPO (ME < 0.049 ppm yr⁻¹ lower, RMSE < 0.305 ppm yr⁻¹, and r>0.888) behaves similarly to the G_ATM derived from the GAW, CTE, and NOAA networks (Figs. 4b and S1). However, noticeable monthly differences between the local and global growth rates, deviating up to approximately 0.8 ppm yr⁻¹, and time shifts are observed (Fig. 4b).

The trend analysis reveals that with development of continental sites, the slope of the trend of annual global CO₂ mole fraction changes from the NOAA network (1.832 ± 0.029 ppm yr⁻¹) to the CTE network (1.859 ± 0.029 ppm yr⁻¹) during 1980–2020 (Fig. S4). However, the G_ATM increased steadily at a rate of 0.030 ± 0.002 ppm yr⁻¹ each year from 1980 to 2020 (Fig. 6a), based on the observations from the three networks (i.e. GAW, CTE, and NOAA). This implies that over long-term periods (here 40 years), the networks with and without continental sites exhibit the same trend in the G_ATM and have little effect on the transient change in the rate of CO₂ increase in the atmosphere. Hence, the role of CO₂ advective transport and mixing in estimating the long-term change in the G_ATM appears negligible. However, a notable difference emerges in the short-term (here 1 month) change in the G_ATM between the networks with and without continental sites (Fig. 6b). El Niño events are known to diminish net global C uptake (due to factors such as droughts, floods, and fires) while increasing global CO₂ growth rate (Sarmiento et al., 2010). During three strong El Niño events, marked as grey bands in Fig. 6b, the G_ATM derived from the GAW and CTE networks (red and green lines) begins to increase before the El Niño events (marked as blue circles in Fig. 6b), approximately 1–2 months earlier than that derived from the NOAA network (black line) and it also reaches its peak during El Niño events (marked as orange circles in Fig. 6b) about 1–2 months earlier (Table S2). This suggests that continental sites can aid in the early detection of G_ATM changes resulting from changes in biogenic emission or uptake. The CTE network (green line) even detects the change 1 month earlier than the GAW network (red line), e.g. for the El Niño 1997–1998 event (Fig. 6b; Table S2). This earlier detection is attributed to the inclusion of even more continental sites in the CTE network (Fig. 1), although the more continental sites also induce the greater variability.

Table 2 presents the global annual CO₂ mole fraction and its G_ATM derived from GAW (GFIT), along with the uncertainty estimates using the bootstrap method. The global average surface CO₂ mole fraction increased from 338.95 ± 0.38 ppm in 1980 to 413.06 ± 0.16 ppm in 2020. Notably, the uncertainty is greater before 1990, primarily due to the limited number of measurement stations worldwide during that period. The average G_ATM for the two decades before 2000 is approximately 1.54 ± 0.08 ppm yr⁻¹. However, in the subsequent two decades, has experienced increases, reaching 1.91 ± 0.05 ppm yr⁻¹ during 2000–2009 and further rising to 2.41 ± 0.06 ppm yr⁻¹ during 2010–2019 (Fig. S5; Table 2).

Table 2 Annual global averaged CO₂ mole fraction (mean in ppm) and its G_ATM (in ppm yr⁻¹) derived from GAW observations using the GFIT method. U(Mean) and U (G_ATM) respectively indicate the uncertainty of the mean and its G_ATM as a 68 % confidence interval. The annual value is averaged over the monthly values of the year.

3.2 Vertical profile of global CO₂ mole fraction

The CTE model simulates CO₂ mole fraction on global 3D grids, enabling us to visualize the modelled vertical CO₂ profile. In the lower atmosphere, highest CO₂ mole fraction is found in the northern mid-latitude region (dark red between 30 and 40^∘ N; Fig. 7a). This area experiences more anthropogenic emissions, which are subsequently transported towards both northern and southern latitudes. The latitudinal and interhemispheric gradient of atmospheric CO₂, as shown in Fig. 7a, is influenced not only by differences in the latitudinal and interhemispheric fossil fuel emissions and seasonal exchanges with terrestrial biota (Denning et al., 1995), but also by atmospheric transport (Patra et al., 2011). As altitude increases, the gradient between the Northern and Southern hemispheres becomes small and levels out at higher altitudes (e.g. >50 km). When comparing the vertical profile change between 2001 and 2020 (Fig. 7b and c), we observe that the CO₂ mole fraction increases slowly in the higher atmosphere (>25 km altitude) compared with the lower atmosphere (<25 km altitude). Figure 7c shows that the vertical gradient (difference between 50 and 0.05 km) changes from approximately 5 ppm in 2001 to around 13 ppm in 2020. The high vertical gradient in 2020 reflects the accumulation of CO₂ in the lower atmosphere, resulting from continuous CO₂ emissions from the surface during 2001–2020 and slow vertical transport. The low vertical gradient in 2001 is partly due to lower surface emissions.

Figure 7 Global vertical profile of CO₂ mole fraction derived from CTE model output. Panel (a) presents the vertical profile in 2020. Panel (b) presents the difference in the vertical profile between 2001 and 2020. Panel (c) presents the annual mean vertical profile from 2001 to 2020. The dots mark CTE vertical level heights and the lines are the linear interpolation between the heights.

Pressure-weighted average CO₂ mole fraction in the lower atmosphere (0–0.35 km altitude) and the entire atmosphere are calculated from CTE output. The annual absolute change in CO₂ mole fraction, computed as the difference between annual means, is more pronounced in the lower atmosphere (orange bars in Fig. S6a) than in the entire atmosphere (blue bars in Fig. S6a). The reason is that the entire atmosphere has a larger air volume than the lower atmosphere, and changes in the surface CO₂ sinks and sources are diluted due to atmospheric horizontal and vertical transport. The CO₂ annual absolute change derived from GAW (GFIT), GAW (WDCGG), and NOAA (represented by red, purple, and brown bars in Fig. S6a) shows small positive or negative differences from the CTE_output (GFIT) and CTE_global (GFIT) across different years. However, over the long term (e.g. on a decadal scale, 2001–2010 and 2011–2020), the CTE model-derived changes in lower and entire atmospheric CO₂ show good agreement (<0.09 ppm yr⁻¹) with the surface observation-based estimate, especially for lower atmospheric CO₂ (<0.07 ppm yr⁻¹). In Fig. S6b, the interannual variability (IAV) of CO₂ mole fraction derived from the CTE model follows a similar temporal pattern as the observation-based IAV derived from the GAW and NOAA network, and especially the IAV of the low-level atmosphere (orange bars) exhibits strong agreement with the observation-based IAV (r>0.971 and RMSE < 0.178 ppm).

3.3 Relationship between the surface CO₂ mole fraction and atmospheric CO₂ mass

The atmospheric CO₂ mass, calculated from the CTE output as a function of air mass and CO₂ concentration (Sect. S3), has increased from 789.46 PgC in 2001 to 877.88 PgC in 2020 (Fig. S7a). The spatial distribution of the atmospheric CO₂ mass is presented in Fig. S7b and c. Monthly global surface CO₂ mole fraction derived from CTE_output (GFIT) and GAW (GFIT), represented as red and blue dots in Fig. 8a, exhibit a similar linear relationship with the monthly atmospheric CO₂ total mass, both showing the same slope of 2.08 ± 0.01 PgC ppm⁻¹. Similarly, NOAA CO₂ (green dots in Fig. 8a) also demonstrates a comparable linear relationship with a slope of 2.09 ± 0.01 PgC ppm⁻¹. Notably, the slopes or conversion factors in Fig. 8a are slightly lower than the factor 2.12 PgC ppm⁻¹ used in Ballantyne et al. (2012) for the period 1980–2010. This minor difference in the conversion factor is expected, considering the different model and data used.

Figure 8 Relationship between the monthly surface CO₂ mole fraction and atmospheric CO₂ mass. The atmospheric CO₂ mass calculated from the 3D CTE output. In (a), the monthly surface CO₂ is derived from the CTE_output (GFIT), GAW (GFIT), and NOAA analysis, presented as blue, red, and green dots, respectively. Panel (b) compares the corresponding interannual variability (IAV) of the atmospheric CO₂ mass and the surface CO₂. The IAV is calculated as the anomaly departure from a quadratic trend.

We further compare the interannual variability (IAV), calculated as the anomaly departure from a quadratic trend, of the atmospheric CO₂ mass and the surface CO₂ (Fig. 8b). The coefficient of the linear relationship closely approaches ∼ 1.0, indicating that the temporal changes in atmospheric CO₂ mass align with the temporal changes in surface CO₂ mole fraction. The CO₂ IAV based on the NOAA network exhibits a slightly closer relationship (r=0.938) with the CTE atmospheric CO₂ mass estimates than the GAW (r=0.861) and CTE (r=0.812) networks. This finding is consistent with the long atmospheric residence time and well-mixed nature of CO₂ in the NOAA network. Overall, the relationship found in Fig. 8 implies that the current surface CO₂ network can effectively serve as an indicator of the CO₂ mass changes throughout the entire atmosphere through a linear relationship.

4 Discussion

Over the past few decades, observational networks have been extended beyond the NOAA MBL network to include more continental sites, such as in the GAW and CTE networks (Fig. 1). These expansions aim to better monitor global CO₂ concentrations and quantify CO₂ sources and sinks. While the continental observations encompass contributions from both substantial sources of anthropogenic emissions and sources/sinks from terrestrial vegetation and soil, these continental observations consistently yield a higher global surface CO₂ mole fraction in the overall global CO₂ analysis, indicating that they are influenced by a bigger net source. We find that the global mean derived from the GAW network is consistently 0.329 (GFIT method) or 0.335 (WDCGG method) ppm higher than that derived from the NOAA network during 1980–2020. Similarly, Tsutsumi et al. (2009) reported a roughly 0.350 ppm higher mole fraction in the GAW network for the years 1983–2006. Notably, the CTE network leads to an even higher global mean (0.422 ppm during 1980–2020), which is likely due to more observational sites located in the Northern Hemisphere, where the highest anthropogenic emissions occur. This also explains the large fluctuation of CO₂ concentrations observed during the winters and summers during 2001–2020 (Fig. 4a). In the future, with the addition of new observation sites, particularly in the Northern Hemisphere, to the existing observational network (e.g. GAW network), we expect that this would lead to higher global surface CO₂ levels and a greater amplitude in the global CO₂ seasonal cycle in the global CO₂ analysis.

Although Friedlingstein et al. (2022) reported a 5.4 % drop (∼ 0.52 PgC) in fossil fuel CO₂ emissions in 2020 (due to restrictions on transport, industry, power, etc., during the COVID-19 pandemic), the increase in annual CO₂ from 2019 to 2020 (2.60 ± 0.16 ppm yr⁻¹) remains at a similar level as from 2018 to 2019 (2.61 ± 0.05 ppm yr⁻¹). In principle, an equivalent drop of roughly 0.25 ppm yr⁻¹ (according to the conversion factor 2.08 PgC ppm⁻¹ in Fig. 8a) or roughly 0.13 ppm yr⁻¹ (according to the annual absolute change; red bars in Fig. S6a) in the growth rate should be visible for the period 2019–2020 due to the declined CO₂ emissions. However, such a short-term human activity induced change in the CO₂ growth rate may be hidden by the natural variability. The bootstrap analysis is used in this study (also in Conway et al., 1994, and Tsutsumi et al., 2009) to estimate the uncertainty of the CO₂ temporal mean and its growth rate and to assess how sensitive the global value is to the distribution of sampling sites. The relatively large uncertainty (± 0.16 ppm yr⁻¹) at the end of 2020 compared with previous years (Table 2) is likely due to an end-effect associated with the curve fitting and filter procedure. The end-effect is a tendency for the growth rate to converge toward the mean value at the end of the record (Conway et al., 1994). Therefore, Conway et al. (1994) suggested that the growth rate curves for the last 6 months should be viewed with caution. Reducing the end effect requires further study, such as using machine learning or bias-correction methods to extrapolate the smoothed trend for a short period (e.g. 1 year) before and after. This extrapolated portion is used exclusively for calculating local mole fraction and growth rate, while it is not included in the global or zonal average, as it could introduce additional uncertainty.

Extrapolation beyond the measurement period extends knowledge gained from a limited period of measurements. During a limited measurement period, we can define the average seasonality, long-term trend, and short-term variation at a measurement site. The long-term trend of an individual site can be extrapolated by various methods, such as referring to the latitude reference time series (Masarie and Tans, 1995) or calculating the mean long-term trend over sites within a certain latitudinal zone (e.g. 30^∘) (Tsutsumi et al., 2009). This extrapolated trend is then combined with the average seasonality to produce estimates beyond the measurement period. However, the extrapolation process relies on the assumption that the relationship of an individual site to the latitude reference remains invariant in time, while in reality the relationship between nearby sites is continuously changing (Masarie and Tans, 1995). In addition, the short-term variation is often ignored or estimated from nearby sites, introducing extra uncertainty into the extrapolation process. In this study, we find that the WDCGG method with extrapolation (GAW (WDCGG+)) results in a global surface CO₂ mole fraction approximately 0.096 ppm higher than the WDCGG method without extrapolation (GAW (WDCGG)) using the same GAW observations, although the extrapolation has a minor effect on the growth rate (Table S1). Therefore, we chose not to use extrapolation beyond the measurement period in our analysis. As the number of long-term measurements increases, the need for such extrapolation becomes less necessary.

Our analysis shows that basing the CO₂ growth rate on GAW surface observations does not introduce a large bias (with an average agreement within 0.016 ppm yr⁻¹) compared with a full atmospheric analysis (Figs. 4b and 5). This full atmosphere CO₂ was provided by the CTE model, in which the global annual mean CO₂ is significantly overestimated compared with GAW observations (e.g. 0.299 ppm higher in CTE_output (GFIT), or 0.186 ppm higher in the CTE_global (GFIT) during 2001–2020). The overestimate derived from the CTE_output (GFIT) is mainly due to more sites in the Northern Hemisphere in the CTE network than in the GAW network. The lower overestimate derived from the CTE_global (GFIT) implies that the biases in CTE outputs are not uniform spatially and tend to balance out. We estimate the CTE bias by comparing the observations and CTE outputs at the same sites, which results in a 0.069 ppm low bias derived from the CTE outputs in calculating the global surface CO₂ mole fraction.

The local growth rate at MLO and SPO generally behaves similarly to the global growth rate derived from the GAW, CTE, and NOAA networks (Figs. 4b and S1). However, the local CO₂ mole fraction and its seasonal cycle noticeably differ from global estimates derived from different observational networks. In this regard, the utilization of individual sites for the evaluation of the global average mole fraction and its growth rate is not precise and can only be used for illustration rather than as a substitute for the proper global average calculation. The local observation sites, often situated away from significant local sources and sinks, such as MLO, provide long-term and high-quality data, serving as reference data for the global CO₂ mole fraction. However, a single observation site cannot capture the CO₂ spatial variability, transport, and mixing. To overcome these limitations, global CO₂ trends and variations are best assessed by integrating data from multiple sources and locations.

Different observational networks (i.e. NOAA, GAW, and CTE) are analysed in this study, revealing differences in calculated global surface CO₂ mole fractions equivalent to the current atmospheric growth rate over a 3-month period. This suggests that the station selection, especially if and how many continental observations are used, has some influence on the derived global surface CO₂ levels, but it is not particularly strong. Nowadays, an increasing number of continental observations are established to monitor biogenic sources and sinks, providing further insight into the climate change and the associated ecosystem processes (Ciais et al., 2005; Ramonet et al., 2020). Such continental observations carry more variability in measurements than the marine observations, which require caution when including them in the mix of stations used to determine global surface CO₂ mole fraction. Our study demonstrates that continental sites can help in early detection of changes in CO₂ growth rate caused by biogenic emission change, such as those resulting from El Niño events. Furthermore, the current observational networks (with and without continental sites) and CTE model show a good agreement on the global CO₂ growth rate, with low or no significant differences within 0.023 ppm yr⁻¹ during 2001–2020 and 0.031 ppm yr⁻¹ during 1980–2020. This implies that the current observation networks (as shown in Fig. 1, representing various ecosystems, sinks, sources, and latitudes) have a similar good capacity to capture changes in the global surface CO₂, although there is the spatial and temporal variability in the CO₂ growth rate (e.g. Conway et al., 1994).

We also notice that the uncertainty in global CO₂ growth rate is approximately 0.07 ppm yr⁻¹, as derived from GAW (GFIT) and averaged over 1980–2020 (Table 2). To reduce the uncertainty to 0.02 ppm yr⁻¹ (equivalent to 1 % of the global CO₂ growth rate), in principle it would theoretically require adding more stations to the current observation network. We conducted an experiment that demonstrates how the uncertainty of the global CO₂ growth rate exponentially increases as the number of land observation sites decreased (Fig. S8). According to our experiment, to achieve the goal of reducing the uncertainty to 0.02 ppm yr⁻¹, 332 land observation sites are required (Fig. S8). However, the required number of sites also depends on their measurement accuracy, consistency, and geographical distribution (i.e. the CO₂ footprint coverage of the observation network and the importance of the network design have been addressed by Storm et al., 2023).

5 Conclusions

The WMO GAW CO₂ network documents the gradual global accumulation of CO₂ in the atmosphere due to human activities. It has been used to assess the large-scale and long-term environmental consequence of fossil CO₂ emission and land use changes. The high-quality observations conducted by the WMO GAW network include not only background stations (most of the NOAA MBL stations) but also continental stations. This comprehensive network enables proper global average calculation. Furthermore, the WMO has initiated a new programme, Global Greenhouse Gas Watch (GGGW), with the aim of establishing a reference network. This network will be built on the high-quality observations already performed under the WMO GAW programme that follows consistent good practices and standards. Although the current monitoring networks have limitations in terms of geographical coverage, data consistency, and long-term measurements, they are well equipped and have the capacity to effectively represent global surface CO₂ mole fraction and its growth rate as well as trends in atmospheric CO₂ mass changes. The three different analysis methods yield very similar global CO₂ increases from 2001 to 2020, which gives us confidence in using any one of them in climate change studies. Continuous monitoring of atmospheric CO₂, based on the current GAW network together with reliable global data integration methods, provides essential information. This includes understanding trends in atmospheric CO₂ concentration, assessing the impacts of past policies, identifying high-emission areas, informing climate models, forecasting future scenarios, and raising public awareness. Policymakers can rely on this information to support their efforts in mitigating global warming.