Enhanced North American carbon uptake associated with El Niño.

Long-term atmospheric CO2 mole fraction and δ13CO2 observations over North America document persistent responses to the El Niño-Southern Oscillation. We estimate these responses corresponded to 0.61 (0.45 to 0.79) PgC year-1 more North American carbon uptake during El Niño than during La Niña between 2007 and 2015, partially offsetting increases of net tropical biosphere-to-atmosphere carbon flux around El Niño. Anomalies in derived North American net ecosystem exchange (NEE) display strong but opposite correlations with surface air temperature between seasons, while their correlation with water availability was more constant throughout the year, such that water availability is the dominant control on annual NEE variability over North America. These results suggest that increased water availability and favorable temperature conditions (warmer spring and cooler summer) caused enhanced carbon uptake over North America near and during El Niño.

INTRODUCTION

Earth’s climate and the carbon cycle are closely linked and mutually affect each other. The El Niño–Southern Oscillation (ENSO) influences weather and climate worldwide (), causing substantial impacts on the rates of global and regional carbon cycling. At the global scale, the interannual variability (IAV) of atmospheric CO2 growth rates shows remarkably strong correlation with ENSO (–), primarily due to ENSO-induced anomalous CO2 fluxes from tropical ocean (, , ) and tropical land (–).

The impact of ENSO on net terrestrial carbon fluxes in extratropical regions was thought to be small, presumably due to weak atmospheric teleconnections between ENSO and extratropical regions and cancellation of compensating effects within subcontinents (, ). Over North America, however, we find from a dense network of atmospheric measurements that atmospheric CO2 mole fractions and δ13CO2 show a broadly consistent response to ENSO (Fig. 1). Here, we detrended and deseasonalized atmospheric CO2 and δ13CO2 data below 1000 m (above ground) at each North American flask air sampling site in the U.S. National Oceanic and Atmospheric Administration (NOAA)’s Global Greenhouse Gas Reference Network (fig. S1 and table S1). We then averaged the residuals to derive average monthly CO2 and δ13CO2 anomalies. The average monthly anomalies of atmospheric CO2 over North America were generally negative around El Niño [Oceanic Niño Index (ONI) ≥ 0.5] and positive around La Niña (ONI ≤ −0.5) periods. They strongly correlated with the ONI with a P-value of 0.0002 (r = −0.38; n = 252). Meanwhile, anomalies of δ13CO2 were anticorrelated with CO2 anomalies (Fig. 1), suggesting that the response of CO2 and δ13CO2 anomalies to ENSO was likely due to changes in terrestrial ecosystem exchange rather than variations in upwind oceanic fluxes. This is because plants (especially C3 plants) take up 12CO2 more efficiently than 13CO2, whereas net oceanic exchange does not strongly discriminate against 13CO2 (versus 12CO2) (, ), so its influence on δ13CO2 is minimal.

Fig. 1

Variability of monthly anomalies of atmospheric CO2 mole fractions and δ13CO2 observations shows a broadly consistent response to ENSO.

(A) The ONI. (B and C) Six-month running averages of monthly CO2 mole fraction and δ13CO2 anomalies averaged across NOAA’s long-term flask air sampling sites over North America (table S1). The number of sites included to calculate the monthly average anomalies of CO2 and δ13CO2 is 7 to 12 for 1995–2003, 16 to 19 for 2004–2007, and 25 to 30 for 2008–2015. Gray shading indicates standard errors of the calculated 6-month running average anomalies. Light yellow shading indicates El Niño periods, whereas light blue indicates La Niña periods.

The variability of atmospheric CO2 mole fraction and δ13CO2 observations over North America reflects an integrated influence of fluxes not only from North America but also from upwind terrestrial ecosystems (i.e., Eurasia). To quantify the extent to which CO2 and δ13CO2 anomalies observed in the atmosphere over North America were associated with changes in North American terrestrial carbon uptake, we inversely modeled North American net ecosystem exchange [NEE = ecosystem respiration (ER) − gross primary production (GPP)] from atmospheric CO2 observations using a high-resolution regional inverse model CarbonTracker-Lagrange (CT-L) (see Materials and Methods) for 2007–2015, a period with extensive ground- and airborne-based air sampling (fig. S1). Along with NOAA’s long-term flask air measurements, we also included additional CO2 observations from NOAA’s continuous tower measurements and measurements made by the Environment and Climate Change Canada and the U.S. Department of Energy available over this period. Estimation of NEE using regional inverse modeling of atmospheric observations requires presubtraction of upwind or background CO2 mole fractions (representative of those in air before it encounters North America) and mole fraction enhancements due to fossil fuel combustion and biomass burning emissions (see Materials and Methods). The remaining spatial and temporal variability in CO2 observations is then used to optimize “a priori” flux estimates such that derived posterior fluxes better represent observed atmospheric mole fractions (than the a priori). In this analysis, we considered three prior flux fields, which were derived from different process-based terrestrial biogeochemical models (TBMs): two Carnegie-Ames-Stanford Approach (CASA) biogeochemical model runs (CASA GFEDv4.1s and CASA GFED-CMS) () and the simple biosphere CASA (SiBCASA) biogeochemical model (fig. S2) ()). We performed an ensemble of 18 inversions (table S2) to include not only uncertainties associated with prior fluxes but also those related to estimated error covariances of prior fluxes and background CO2 mole fractions (see Materials and Methods and table S2) that cannot be explicitly derived from a single inversion.

RESULTS AND DISCUSSION

IAV of North American NEE

The mean annual net North American land CO2 flux derived from CT-L (posterior NEE + prescribed emissions from fires) was −0.70 (−1.24 to −0.34) PgC year−1 between 2007 and 2015, considering the range of 18 inversion results plus the 2σ uncertainties from each inversion. This agrees well with estimates from the global inverse models such as CarbonTraker (CT2016) (−0.56 PgC year−1) and CarbonTracker-Europe (CTE2016) (−0.67 PgC year−1) over these same years. It is also within the range of earlier estimates provided by other inverse models between 2000 and 2015 (fig. S3) (–).

Monthly NEE anomalies are used to diagnose IAV and were estimated from each inversion by subtracting the derived monthly NEE from the monthly fluxes averaged for the same months between 2007 and 2015. Monthly NEE anomalies discussed below are average anomalies across 18 inversions. Derived monthly NEE anomalies varied by as much ±1.5 PgC year−1 during 2007–2015 and are dominated by flux anomalies derived for temperate North America (~ 25° to 50°N) (Fig. 2A). The temporal variability in derived monthly NEE anomalies is similar to the variability in measured atmospheric CO2 mole fraction and δ13CO2 (Fig. 1) and anticorrelated with the ONI with 1 to 3 months lag (r < −0.36; P < 0.03). The temporal variability and spatial patterns of the derived NEE anomalies suggest broad responses of both temperate and boreal North American terrestrial ecosystems to ENSO (Fig. 2 and fig. S4) with larger anomalies being derived for three biomes: forest field, crops, and coniferous forests (fig. S4). These forested biomes have previously been identified as the primary contributors to the overall IAV of North American NEE ().

Fig. 2

Anomalies of derived NEE over North America display distinct responses to different phases of ENSO between 2007 and 2015.

(A) Derived monthly NEE anomalies from 18 inversions considered in this study from CT-L (gray shading), their monthly means (black thin line), and 6-month running averages (black thick line). Blue and red solid lines indicate the 6-month running averages of monthly NEE anomalies derived for boreal and temperate North America. Light yellow shading indicates El Niño, whereas light blue indicates La Niña. (B) Average monthly NEE anomalies during El Niño, neutral, and La Niña periods [which were indicated by yellow, white, and blue shadings in (A)], derived from TBMs (CASA GFEDv4.1s, CASA GFED-CMS, and SiBCASA), CarbonTracker (CT2016), and CT-L (this study). Anomalies derived from CT-L were indicated by the boxplot with red lines representing the medians, blue boxes indicating the 25th and 75th percentiles, and black dashed bars indicating the minimums and maximums of the ranges calculated from the 18 inversions.

Our derived NEE anomalies from CT-L show that, during El Niño, there was an additional 0.61 (0.45 to 0.79) PgC year−1 (the range from 18 inversions) sequestered by North American terrestrial ecosystems compared to the La Niña periods (Fig. 2B). This response seems to be asymmetric, with a stronger response to the El Niño phase [−0.38 (−0.48 to −0.30) PgC year−1] than to the La Niña phase [0.22 (0.07 to 0.31) PgC year−1] (Fig. 2B). The natural variability of North American NEE resulting from climate variability largely associated with ENSO is much larger than the IAV of fossil fuel or fire emissions (1σ ≤ 0.05 PgC year−1 between 2007 and 2015). It is also substantial compared to the magnitude of natural or anthropogenic fluxes of greenhouse gases over North America: It is equivalent to roughly 80% of North American annual NEE (−0.78 PgC year−1), one third of annual total North American fossil fuel CO2 emissions (1.7 PgC year−1) (, ) and twice as large as the total U.S. anthropogenic non-CO2 greenhouse gas emissions (CH4, N2O, and fluorinated gases) (0.3 PgC-equivalent year−1) (, ). This implies that shifts in climate (i.e., ENSO phases) have the potential to largely alter the overall carbon fluxes from North America.

The response of North American net carbon uptake to ENSO identified here is not present in the a priori fluxes that were calculated from process-based TBMs (fig. S5). It only emerges in posterior fluxes that were optimized using atmospheric observations (fig. S5). Fluxes estimated from the global inverse model CarbonTracker (CT2016) (Fig. 2B and fig. S5) (), which used a similar suite of atmospheric CO2 observations and the same prior fluxes (CASA GFEDv4.1s and CASA GFED-CMS) as CT-L, also display a clear response to ENSO but with a magnitude less than half of that estimated by CT-L (Fig. 2B and fig. S5). The difference in North American NEE response to ENSO derived from CT-L and CT2016 may arise from the higher resolution of CT-L in both the atmospheric transport simulation and inverse modeling setup. In CT2016, air transport was simulated at 1° × 1° resolution over North America (between 20° to 64°N and 132° to 60°W and at 2° latitude × 3° longitude over the rest of the globe), and scaling factors of fluxes were derived on a weekly basis for each ecoregion. In CT-L, we computed air transport at 10-km resolution over most of North America (see Materials and Methods) and optimized weekly scaling factors of fluxes at 1° × 1° resolution with adjustment of the average weekly diurnal cycle of NEE (see Materials and Methods). Thus, compared to CT2016, CT-L has more flexibility to adjust diurnal cycles and spatial patterns of NEE. High-resolution (1° × 1° × 3 hourly) NEE estimates derived using the same transport simulations as CT-L from 2007 to 2012 but using a geostatistical inverse modeling approach () show an almost identical magnitude of response to ENSO over this period as our estimate (fig. S5B). This supports the hypothesis that the difference in the estimated NEE response to ENSO from CT-L and CT2016 may result from the difference in transport models and inverse modeling resolution. Compared to the global CarbonTracker, fluxes derived from our regional model are likely more accurate, as indicated by better agreement with the atmospheric CO2 observations at a 99.9% confidence interval (figs. S6 and S7). This may be because the use of higher-resolution models improves air transport simulation (, ) and reduces spatial and temporal aggregation errors when deriving fluxes (figs. S8 and S9) (, ). The response of North American carbon uptake to ENSO likely also occurred before the availability of our flux estimates that begin in 2007, as indicated by atmospheric CO2 mole fraction and δ13CO2 observations (Fig. 1) and in fluxes derived from CarbonTracker (fig. S10).

Climate Drivers for IAV of North American NEE

To understand why North American NEE exhibits such a consistent response to ENSO, we first investigated how climate variables drive the IAV of North American net carbon uptake. Current understanding of these relationships is so poor that projections of future carbon sequestration rates disagree in both sign and magnitude (). In particular, the relative roles of air temperature and water availability in controlling NEE have been intensely debated [e.g., (, , –)]. Jung et al. () partially resolved this debate by suggesting that spatially compensatory effects of water-dominated controls on local scales minimize its influence so that air temperature is the dominant control on global yearly land carbon sink. Over temperate North America, precipitation was recently suggested to dominate control of GPP and total ER (). However, because of their compensating effects, it remains unclear what dominates the control on the IAV of NEE. Furthermore, studies have also suggested that the control of climate variables on IAV of North American NEE may vary geographically (, ) and seasonally (). Therefore, the relative importance of air temperature versus water availability on controlling North American carbon uptake needs to be further explored.

With our newly derived NEE from atmospheric CO2 observations, we conducted a correlation analysis between monthly NEE anomalies and anomalies of climate variables, including air temperature, precipitation, relative humidity (RH), water vapor pressure deficit (VPD; deficit between saturation water vapor pressure and actual water vapor pressure in ambient air), and soil moisture (SM). We find that North American NEE anomalies correlated significantly with anomalies of hydrological parameters such as VPD, RH, SM, and precipitation over temperate North America at a 95% confidence interval (Table 1 and table S3). Although the correlation with precipitation is statistically significant, it is much weaker than with other hydrological parameters on the continental scale. In addition, NEE anomalies are significantly correlated with anomalies of VPD, RH, SM, and precipitation only when a 1- to 2-month lag was considered in NEE (table S3), suggesting a time lag for the response of net carbon fluxes to hydrological conditions. This has also been noted over the Amazonian rainforest ().

Table 1

Correlation coefficients between CO2 flux anomalies over North America and anomalies of area-weighted average air temperature, precipitation, RH, VPD, and SM over temperate North America, using prior fluxes (CASA GFEDv4.1s) and derived posterior fluxes from CT-L.

Correlations were calculated with fluxes lagging climate variables by 1 month. The P value associated with each correlation is included in parentheses. The correlation is higher for yearly average anomalies than that for monthly anomalies as noise in the datasets are smoothed out.

Climate variables	Prior fluxes (monthly)	Posterior fluxes (monthly)	Posterior fluxes (yearly)
Air temperature	−0.09 (0.34)	−0.01 (0.90)	0.30 (0.43)
Precipitation	−0.03 (0.78)	−0.21 (0.03)	−0.78 (0.01)
RH	−0.14 (0.15)	−0.38 (<0.0001)	−0.69 (0.04)
VPD	0.16 (0.09)	0.39 (<0.0001)	0.75 (0.02)
SM	−0.23 (0.02)	−0.32 (0.001)	−0.66 (0.05)

When considering all months of the year together, monthly North American NEE anomalies are not correlated with air temperature anomalies (Table 1 and table S3). However, when considered on a seasonal basis, correlations emerge that change signs with season (Fig. 3): In spring (March to May), a negative correlation exists between temperature and NEE anomalies, suggesting that warm conditions in spring enhance North American carbon uptake; whereas in summer (June to August), the opposite correlation exists, implying that increased summertime temperatures reduce terrestrial carbon uptake. During fall (September to November) and winter (December to February), the correlation between air temperature and NEE anomalies is statistically insignificant. Note that the spring and summer correlations were significant with or without a 1-month lag. These correlations are even stronger than those between hydrological variables and NEE anomalies in these seasons (Fig. 3). However, in contrast to opposite correlations between seasons (in the relationship between air temperature and NEE anomalies), correlations between hydrological variables and NEE anomalies were relatively more constant throughout the year (Fig. 3).

Fig. 3

Climate drivers for anomalies of North American NEE and their response to ENSO.

(A) Correlation coefficients between NEE anomalies and anomalies of air temperature, precipitation, VPD, and SM in different seasons using prior fluxes (computed from TBMs) and fluxes derived from CT-L. Strong correlations with P < 0.1 are indicated by color-filled symbols, whereas weaker correlations with P > 0.1 are shown in empty symbols. (B) Differences of average air temperature anomalies between El Niño and non–El Niño (La Niña and neutral conditions) periods in spring and summer and difference of average VPD anomalies between El Niño and non–El Niño periods averaged for all seasons. While the impact of El Niño on North American air temperature is opposite between spring and summer, its impact on North American VPD (and other hydrological variables) is relatively more constant throughout the year (fig. S14).

The correlations between NEE anomalies and hydrological variables discussed above were mostly insignificant in the prior fluxes computed from different process-based TBMs, CASA (CASA GFEDv4.1s and CASA GFED-CMS), and SiBCASA (Fig. 3, Table 1, and table S3). The correlations emerge only from the optimization of these fluxes through inverse modeling of the atmospheric CO2 observations (Fig. 3, Table 1, and table S3). This suggests an underestimation of the sensitivity of NEE to water availability or moisture conditions over North America by TBMs. With no time lag, the opposite seasonal temperature effects of air temperature on NEE anomalies in spring and summer derived in the posterior fluxes of CT-L did not exist in prior TBM-computed fluxes (Fig. 3). When considering a 1-month lag, anomalies of prior fluxes also showed an opposite correlation with temperature anomalies, but much weaker than derived from the posterior fluxes. Poor representation of the carbon uptake and climate relationships by the TBMs likely partially explains why the computed North American fluxes from TBMs have not shown a response to ENSO. The poor representations of IAV of North American carbon uptake and their relationships with climate variables in the TBMs have also been noted elsewhere (, , ).

Correlations between NEE anomalies and climate variables over the major ecoregions (that have larger net carbon uptake; fig. S4) are similar to those noted on the continental scale (fig. S11). They show a strong and opposite temperature control on NEE in spring and summer, although hydrological conditions (VPD, SM, RH, and precipitation) were also important especially during summer in controlling the IAV of NEE. During fall, hydrological variables seem to be more important in controlling net carbon uptake than temperature in many ecoregions, except over some biomes in boreal North America: semi-tundra and northern taiga (fig. S11), where air temperature positively correlated with NEE anomalies. This relationship was noted previously and suggests net carbon losses from northern high-latitude ecosystems in response to autumn warming (). In winter, the correlations between NEE and climate anomalies are the weakest among all seasons with precipitation and moisture conditions more important in controlling ecoregional NEE than air temperature.

Spatial distribution of Pearson’s correlation coefficients between anomalies of hydrological variables and NEE indicates stronger control of water availability on NEE over temperate than boreal regions (fig. S12). The consistent sign in the Pearson’s correlation coefficients across temperate North America indicates small spatial cancellation effects over this region, as opposed to strong spatial cancellation noted on a global scale (). For the correlation between air temperature and NEE anomalies (fig. S13), it was stronger in boreal and eastern temperate North America than that in western temperate North America in spring, whereas in summer, they were most strongly correlated across the temperate region, although the correlations between summer and spring had opposite signs. Climate sensitivities, as derived from slopes of a linear regression between climate variables and NEE anomalies, suggest generally larger sensitivities over temperate than boreal regions (figs. S12 and S13), likely due to much larger variability of NEE in temperate ecosystems (Fig. 2).

UNDERSTANDING THE RESPONSE OF NORTH AMERICAN NEE TO ENSO

Earlier studies have shown that, around El Niño periods, there was above normal precipitation in temperate North America from 1875 to 1990 (, ). We found that a similar impact was also evident for the last several decades (figs. S14 and S15). The strength of ENSO, expressed by the ONI, not only strongly correlates with anomalies of precipitation but also significantly correlates with anomalies of RH, VPD, and SM over temperate North America (with ± a few-month time lag) (table S4 and fig. S15), suggesting more precipitation, higher RH, smaller VPD, and greater SM over temperate North America before, during, and after an El Niño event (figs. S14 and S15). This ENSO-climate correlation over temperate North America appears to be much stronger after 1980 than during 1950–1980 (fig. S16).

Besides increased water availability (increased precipitation, SM, and RH and smaller VPD), favorable temperature conditions further enhanced North American carbon uptake around El Niño over the last decade. In summer, it was cooler during El Niño than non–El Niño (neutral and La Niña) periods (Fig. 3 and fig. S14). As temperature positively correlates with NEE in summer, decreased air temperature would result in more negative NEE (more carbon uptake). Meanwhile, warmer than usual spring was observed during El Niño, especially over the boreal region (Fig. 3 and fig. S16), and would further enhance terrestrial carbon uptake over North America, as air temperature negatively correlates with NEE anomalies (Fig. 3).

The response of North American terrestrial ecosystems to ENSO derived here is opposite to the ENSO–carbon flux relationship observed over tropical terrestrial ecosystems (, , ). In the tropics, the climatic impact of El Niño tends to be reduced precipitation and increased air temperatures (drought) (), whereas the imprint of El Niño over North America appears to be wetter with increased temperatures in winter and spring and decreased temperatures in summer (fig. S14). When compared to the increased net biospheric carbon flux over the tropics estimated for the recent ENSO event (2.5 ± 0.3 PgC year−1 increase between the May 2015 to April 2016 El Niño period and the 2011 La Niña period) (), we estimate an offsetting carbon uptake in the extratropics from North America that is roughly a quarter [24% (16 to 36%)] as large, considering the anomalies we estimated for 2007–2015 [−0.61 (−0.79 to −0.45) PgC year−1]. Our approximation neglects the difference in time periods associated with both estimates. Neither did we consider the ENSO response of other extratropical land regions (e.g., Eurasia), which may additionally enhance this offset. Our results highlight the importance of improving the understanding of regional carbon-climate relationships, which represent a major uncertainty in future climate projections (). Furthermore, since climate change will likely result in more frequent and intense El Niño events in the future (), our results also underscore the importance of improving quantification of regional response of carbon cycling rates to climate anomalies associated with ENSO to more accurately estimate the overall climate impact of future ENSO events.

MATERIALS AND METHODS

Atmospheric observations

Atmospheric CO2 observations used in this analysis for 2007–2015 were from the GLOBALVIEWplus v2.1 ObsPack (www.esrl.noaa.gov/gmd/ccgg/obspack/) prepared by the Global Monitoring Division (GMD), Earth System Research Laboratory (ESRL), NOAA, United States. Those observations included flask air measurements (table S1) and continuous in situ measurements primarily made by GMD with additional contributions from the Environment and Climate Change Canada and U.S. Department of Energy (fig. S1). Flask air samples were collected from near the Earth’s surface, towers, and aircraft, whereas continuous measurements were made only from towers. CO2 observations included in inversions between 2007 and 2015 were ground-based observations selected and simulated by the NOAA CarbonTracker (CT2016), plus aircraft data that had above zero sensitivity to upwind surface fluxes (or “footprints”) (). The selection criteria for atmospheric CO2 observations used in CarbonTracker’s assimilation are described on the NOAA CarbonTracker website (www.esrl.noaa.gov/gmd/ccgg/carbontracker/CT2016_doc.php#tth_sEc7.1). The total number of observations used in the inversions conducted by this study was 88949, 16.9% of which were from aircraft. CO2 observations used for analyses before 2007 were only based on measurements made from flask air samples collected by GMD.

Atmospheric δ13CO2 observations were made from flask air samples collected by GMD (fig. S1 and table S1) and analyzed by the Institute of Arctic and Alpine Research (INSTAAR), University of Colorado. Methods for measurements of δ13CO2 can be found in White et al. () and Trolier et al. (). INSTAAR used a local realization of the VPDB-CO2 scale based on many measurements made in the early 1990s (). Since then, the local scale has been propagated by CO2-in-air cylinders. Although there were documented offsets from other laboratories, long-term surveillance cylinders showed the consistency of the INSTAAR scale over time.

Anomalies of atmospheric CO2 and δ13CO2 were calculated from measured atmospheric CO2 dry air mole fractions or δ13CO2 subtracted from their long-term trends and multiyear averaged seasonal cycle at each site, computed from the algorithms given by Thoning et al. (). The detrended and deseasonalized residuals were then averaged within each month to obtain monthly average anomalies at individual sites.

Estimating North American NEE using a high-resolution regional inverse model: CT-L

We recently developed a high-resolution regional inverse modeling framework for estimating North American regional fluxes by assimilating atmospheric data: CT-L (www.esrl.noaa.gov/gmd/ccgg/carbontracker-lagrange/). CT-L used the high-resolution Weather Research and Forecasting–Stochastic Time-Inverted Lagrangian Transport (WRF-STILT) model () to simulate atmospheric transport and compute footprints. The WRF model fields had 10-km spatial resolution below 55°N over temperate North America and 30- to 40-km spatial resolution above 55°N during 2007–2010. After 2010, the high-resolution (10 km) domain was extended to 55°N, covering most of boreal North America except for Alaska. STILT footprints have 1° × 1° × 1 hour resolution and represent simulated upwind influences over 10 days before each measurement. With precomputed footprints, we estimated NEE using Bayesian inversion algorithms described below and in our earlier papers (, , ). NEE here is defined as the difference between ER and GPP. Negative NEE indicates CO2 uptake by terrestrial biomes from the atmosphere, and positive NEE implies CO2 releases from terrestrial biomes to the atmosphere.

We assumed that each CO2 mole fraction observation (zCO) can be modeled as the sum of the CO2 mole fraction in background air upwind of the North American continent (zbg) and recent contributions from NEE (zbio), fossil fuel emissions (zff), and emissions from wildfires (zfire). Background CO2 mole fraction (zbg) for each observation was estimated by three approaches using CarbonTracker (CT2016) simulated four-dimensional (4D) CO2 mole fraction field, an alternative 4D CO2 mole fraction field estimated from an empirical method, and back trajectories (see detailed description below). Contributions from fossil fuel (zff) and fire (zfire) emissions were estimated by multiplying footprints computed from WRF-STILT with fossil fuel and fire emission inventories. We averaged the two fossil fuel emission products (“Miller” and “ODIAC” datasets) and fire emission products (“GFED4.1s” and “GFED_CMS”) used in CarbonTracker (CT2016) on a 3-hourly basis and used the average 1° × 1° × 3 hour emissions to compute their associated CO2 mole fraction enhancement. The remaining variability of CO2 observations (zbio) (the observed CO2 mole fractions minus estimated background mole fractions and mole fraction enhancements due to fossil fuel and fire emissions) was then used in our Bayesian inversion framework to estimate North American NEE.

The cost function for Bayesian inverse modeling can be written aswhere zbio = zCO − zbg − zff − zfire; K is the estimated change in the CO2 mole fraction due to fluxes over 10 days before each measurement using prior 3-hourly NEE (Pbio) or K = HPbio; H is the sensitivity of atmospheric observations to upwind fluxes or footprints as mentioned above. λ is a vector of prior scaling factors with values of ones. λ is a vector of optimized posterior scaling factors for prior fluxes that were solved through inversions. We solved for eight scaling factors (λ) per week for each 1° × 1° grid cell, one for each 3-hourly time of day (i.e., a weekly average adjustment to the diurnal cycle). Prior NEE (Pbio) had 1° × 1° × 3-hour resolution, so posterior NEE (Fbio; Fbio = λPbio) did account for day-to-day variability within a given week to the extent that the prior fluxes captured variability such as ecosystem responses to local weather conditions. The length of λ or λ was equal to the number of fluxes being estimated, which equated to the number of grid cells multiplied by the number of weeks in each batch inversion that spanned a year plus 2 weeks before and after the year. R and Q are covariance matrices of model-data mismatch errors and prior relative flux errors, both of which were estimated by maximum likelihood estimation (MLE) ().

The posterior scaling factors of the 1° × 1° NEE fluxes were solved according to Eq. 2, while the posterior error covariance matrix (C) was computed using Eq. 3 (A and B):where V is the posterior error covariance matrix of the 1° × 1° scaling factors. The 1σ error of the 1° × 1° posterior fluxes is equal to the square root of the diagonal elements of C. Note that it is impractical to compute the V and C matrices explicitly due to their large size, so we adapted the methods described by Yadav and Michalak () to calculate monthly or annual V and C. The 1σ errors of posterior fluxes aggregated to ecoregions () or for North America, temperate or boreal North America (defined in fig. S1), are equal to the square root of the sum over rows and columns of C corresponding to the grid cells representing those regions, considering their cross correlations. The definition of ecoregions and temperate or boreal North America used in this analysis is consistent with that used in CarbonTracker (). The errors mentioned above are statistical errors computed from each inversion run. Additional uncertainties associated with prior fluxes, estimated error covariance matrices, and estimate background mole fractions were considered by performing an ensemble of 18 inversions with different choices for those parameterizations (table S2). Systematic errors in air transport simulation were not specifically examined in this study. However, the IAV, which we focused on in this paper, is expected to be relatively insensitive to systematic errors of air transport as found by previous studies (, , ). Similarly, errors in fossil fuel combustion and biomass burning (fire) emissions may be propagated into the estimate of NEE. Because the IAV of fossil fuel combustion and biomass burning (fire) emissions are relatively small, the contribution of their associated errors to the estimated IAV of NEE is likely negligible.

Prior net ecosystem exchange

We used three different prior NEE (fig. S2) to test the sensitivity of derived posterior fluxes to prior flux magnitudes, amplitudes of diurnal variability, and spatial distributions. The three prior fluxes are: (i) and (ii) two prior flux products used in CT2016, which estimated 3-hourly NEE fluxes from monthly NEE calculated from two different CASA biogeochemical model runs (CASA GFED-CMS and CASA GFEDv4.1s) () and (iii) fluxes computed from the SiBCASA terrestrial carbon cycle model (). The diurnal cycle used for the CT2016 prior fluxes was the result of temporal downscaling according to algorithms developed by Olsen and Randerson (), whereas in SiBCASA, the diurnal flux was simulated explicitly. The SiBCASA 3-hourly prior fluxes have a diurnal cycle with amplitude approximately half as large as CT2016 prior fluxes (CASA GFED-CMS and CASA GFEDv4.1s) (fig. S2). There were also differences in spatial distributions and magnitudes of monthly and annual NEE between CT2016 prior fluxes and SiBCASA fluxes.

Estimating error covariances

Error covariances (R and Q) are critical parameters that determine the weight between prior fluxes and atmospheric observations in the final posterior solution. Many previous studies have relied on “expert judgment” to specify error covariances. Here, we used a more objective approach, MLE (, ), for estimating R and Q. This approach guarantees a reduced χ2 of 1 in the posterior solution. In this study, we estimated a seasonally (every 3 months) and interannually varying model-data mismatch error (the square root of the diagonal element in R) on a site-by-site basis and a seasonally and interannually varying prior flux error on a relative scale (the square root of the diagonal element in Q) across North American land. The MLE-estimated model-data mismatch errors are typically in a range of 0.4 to 3.4 ppm (0 to 90th percentile) and estimated prior flux errors are ~100 to 300% for the three priors used in inversions, CASA GFED-CMS, CASA GFEDv4.1s, and SiBCASA, depending on seasons or years. Since measurement errors are typically on an order of 0.1 ppm (), the MLE-estimated model-data mismatch errors are likely dominated by the other components, such as transport errors, model representation errors, and background errors.

We also used MLE to optimize the e-folding correlation scales (in time and space) in each batch inversion for each year. Here, we assumed the correlation decays exponentially in time and space as applied previously in other inversion analyses [e.g. (, )]. The estimated e-folding correlation scales in prior flux errors are in an approximate range of 20 to 50 days and 200 to 600 km depending on the year or the prior flux (table S2). To test the sensitivity of modeled posterior fluxes to the prior flux error correlation scales, we also performed inversions with assumed e-folding correlation scales (7 days and 1000 km) that were much different from what were derived from the MLE method. Note that, for the temporal correlation scales, we applied correlations across days at the same time of day (e.g., 1:00 p.m. today is correlated with 1:00 p.m. tomorrow) and no within-day correlation (e.g., 1:00 p.m. today is not correlated with 3:00 p.m. or 11 a.m. today). This approach was also used by Gourdji et al. (). As expected, our inversions using longer spatial correlation scales (1000 km versus 200 to 600 km) resulted in spatially smoother fluxes, and inversions with shorter temporal correlation scales (7 days versus 20 to 50 days) resulted in larger temporal variability in derived monthly fluxes.

Background CO2 mole fraction estimation

Background CO2 mole fraction for each CO2 observation was estimated using three approaches combining 4D background CO2 mole fraction fields with WRF-STILT back trajectories: (i) The approach “CTBG” used CT2016’s simulated 4D CO2 mole fraction field and is the average value of the sampled CO2 mole fractions from this field at locations where the 500 back trajectory particles exited the WRF modeling domain. Any trajectory particles remaining in the domain after 10 days were sampled at their ending locations within the domain. (ii) The approach “EBG” used an alternative 4D CO2 mole fraction field that was constructed from observations made in the remote atmosphere either within the marine boundary layer at remote sites or within the free troposphere at aircraft sites over North America. The free troposphere reference surfaces were created for altitudes > 3 km above sea level, using forward and backward trajectory ensemble to extend the influence of individual aircraft observations throughout the domain. An extended aircraft dataset was created by assuming the measured value was valid along these trajectories at 3-hourly intervals. Trajectories were terminated when any trajectory particle dropped below 3 km above ground level, as surface fluxes would be expected to modify the value at low altitudes. The extended dataset has the highest density in the vicinity of the actual measurement location, and the density decreases as the trajectories disperse. The extended dataset was gridded to 5° × 5° × 1 km (longitude × latitude × altitude intervals), and smooth curves were fitted to each bin using methods described by Thoning et al. (). Random errors were estimated on the basis of the curve fit residuals and were scaled according to the number of actual observations (instead of the extended observations). The free troposphere reference surfaces were combined with separate Atlantic and Pacific marine boundary layer reference surfaces (constructed using the same methodology used to construct the global marine boundary layer reference, www.esrl.noaa.gov/gmd/ccgg/globalview/). The free troposphere reference surfaces and Atlantic and Pacific marine boundary layer reference surfaces were combined to provide time-dependent lateral and upper boundary values for North America. WRF-STILT trajectories for each observation used in the inversions were sampled at their endpoints as for the method (i). Any trajectories terminating within the continental boundary layer were assigned a latitude-dependent Pacific or Atlantic marine boundary layer value depending on whether the trajectory endpoint longitude was west or east of 100°W. The EBG approach has the advantage of being directly constrained by aircraft and marine boundary layer observations, but it does not capture synoptic variability and it lacks an accurate representation of the continental boundary layer. For midcontinent and eastern sites, up to 20% of trajectories terminated within the continental boundary layer so biases may be nonnegligible. (iii) Considering known seasonal biases in CT2016’s CO2 simulations compared to observations (fig. S6) and the limitations of the EBG approach, the third approach applied a correction for background mole fractions estimated from CTBG for each season each year on a site-by-site basis (CTBGcorr). For this correction, we first selected a subset of “background observations” that had minimal sensitivity [summed footprints ≤ 1 ppm (μmol m−2 s−1)−1] to land fluxes over the North American land (fig. S17). Then, a correction was made on the basis of an average seasonal difference between CTBG and the selected background observations below 3 km above ground at each site (fig. S18). This approach was applied previously to anthropogenic gases ().

Observing system simulation experiments using pseudo-observations to optimize inversion configuration

We conducted a suite of observing system simulation experiments (OSSEs) to determine an optimal inversion configuration that yields reliable flux estimation with a reasonable computing efficiency. In these experiments, we first assumed 3-hourly SiBCASA fluxes were the “truth” and computed pseudo-observations with added random noise (1σ = 0.1 ppm). The computed pseudo-observations had the same sampling times and locations as real observations. We then used one of CT2016’s prior 3-hourly fluxes (CASA GFEDv4.1s) as the prior of these experiments. We computed the scaling factors of prior fluxes at three different spatial and temporal resolutions using Bayesian inversion algorithms described above: (i) Similar to the inversion setup of the CarbonTracker, we solved for weekly scaling factors of prior fluxes on ecoregional scales [ecoregions are defined in the same way as CarbonTracker () and shown in fig. S1]. (ii) Compared to the configuration (i), we increased the spatial resolution of the inversion setup by solving for weekly scaling factors for fluxes on 1° × 1° grid scales. (iii) Compared to configuration (i), we increased both spatial and temporal resolution of the inversion by solving for weekly scaling factors at 1° × 1° with an additional optimization on weekly average diurnal cycles on 1° × 1° grid scales (eight scaling factors per week per grid cell). The third setup is the one that was used in real-data inversions in CT-L. To rule out the impact of different covariance parameters on the inversed results, we applied the same covariance parameters calculated from MLE to all three configurations. Note that these OSSEs represent idealized cases with no systematic biases and no transport errors considered. However, they are useful for testing the mechanics of the inversion and optimizing the mathematical construct of the inversion. If an inversion setup does not work to retrieve true fluxes in the OSSEs, then it would be even harder to derive accurate fluxes in real-data inversions.

When we solved for weekly scaling factors of prior fluxes on 1° × 1° grid scale with an adjustment on the diurnal cycle of prior fluxes, the derived posterior fluxes had the smallest aggregation errors in time and space, compared to those obtained from weekly scaling on ecoregional or grid-scale fluxes described in configurations i and ii above (figs. S9 and S10). They also better agreed with true fluxes on ecoregional scales, especially over regions such as semi-tundra and northern taiga. Allowing for adjustment of diurnal cycles of NEE has already been implemented in previous regional CO2 inversion studies (, , , ) but has typically been neglected in global CO2 inversions.

Climatological data

The ONI between 1950 and 2017 was downloaded from NOAA’s National Weather Service Climate Prediction Center (http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). It represented 3-month running mean of ERSST.v5 sea surface temperature anomalies in the Niño 3.4 region (5°N to 5°S, 120° to 170°W). An El Niño condition was defined as ONI ≥ 0.5. All periods with ONI < 0.5 were considered as non–El Niño conditions, including neutral (−0.5 < ONI < 0.5) and La Niña (ONI ≤ −0.5) conditions.

Gridded monthly air temperature, precipitation, and RH were from NOAA’s National Centers for Environmental Prediction’s North American Regional Reanalysis (NARR), provided by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website at www.esrl.noaa.gov/psd/. NARR datasets were available at 32-km resolution over North America for 1979 to 2017. Gridded (0.5° × 0.5°) monthly SM data were from the NOAA Climate Prediction Center’s SM database. They were downloaded from the NOAA/OAR/ESRL PSD (www.esrl.noaa.gov/psd/data/gridded/).

VPD was calculated from saturation water vapor pressure (Vsat) at ambient air temperature and the RH in the ambient airwhere, Vsat was computed fromT is the air temperature in degrees Celsius and Vsat is in kilopascal.

Monthly average air temperature, precipitation, RH, SM, and VPD over boreal and temperate North America were calculated on the basis of area-weighted averages over grid cells in those regions. To examine the monthly anomalies that resulted from IAV, we first removed the decadal and multidecadal variabilities by subtracting their 10-year running means from the monthly average data. We then removed the average seasonality from the monthly datasets to derive monthly anomalies of air temperature, precipitation, RH, SM, and VPD.

Correlation between the environmental variables (air temperature, precipitation, RH, SM, and VPD) and ONI was performed on the basis of their 3-month averages. Since ONI data were 3-month running averaged data, we only extracted values that did not have a temporal overlap; for each year, we only had 4 independent ONI values instead of 12 ONI values. This approach allowed us to correctly count the degrees of freedom in the climate datasets and prevented overestimating their correlations.