Increasing summer net CO 2 uptake in high northern ecosystems inferred from atmospheric inversions and comparisons to remote-sensing NDVI

Warmer temperatures and elevated atmospheric CO2 concentrations over the last several decades have been credited with increasing vegetation activity and photosynthetic uptake of CO2 from the atmosphere in the high northern latitude ecosystems: the boreal forest and arctic tundra. At the same time, soils in the region have been warming, permafrost is melting, fire frequency and severity are increasing, and some regions of the boreal forest are showing signs of stress due to drought or insect disturbance. The recent trends in net carbon balance of these ecosystems, across heterogeneous disturbance patterns, and the future implications of these changes are unclear. Here, we examine CO2 fluxes from northern boreal and tundra regions from 1985 to 2012, estimated from two atmospheric inversions (RIGC and Jena). Both used measured atmospheric CO2 concentrations and wind fields from interannually variable climate reanalysis. In the arctic zone, the latitude region above 60 N excluding Europe (10W–63 E), neither inversion finds a significant long-term trend in annual CO2 balance. The boreal zone, the latitude region from approximately 50–60 N, again excluding Europe, showed a trend of 8–11 Tg C yr−2 over the common period of validity from 1986 to 2006, resulting in an annual CO2 sink in 2006 that was 170–230 Tg C yr−1 larger than in 1986. This trend appears to continue through 2012 in the Jena inversion as well. In both latitudinal zones, the seasonal amplitude of monthly CO2 fluxes increased due to increased uptake in summer, and in the arctic zone also due to increased fall CO2 release. These findings suggest that the boreal zone has been maintaining and likely increasing CO2 sink strength over this period, despite browning trends in some regions and changes in fire frequency and land use. Meanwhile, the arctic zone shows that increased summer CO2 uptake, consistent with strong greening trends, is offset by increased fall CO2 release, resulting in a net neutral trend in annual fluxes. The inversion fluxes from the arctic and boreal zones covering the permafrost regions showed no indication of a large-scale positive climate–carbon feedback caused by warming temperatures on high northern latitude terrestrial CO2 fluxes from 1985 to 2012.


Introduction
The high northern latitudes, including the tundra and boreal forest regions, are particularly vulnerable to the effects of climate change as this region has been experiencing dramatic changes in recent climate.Warming in northern ecosystems results in many physical and ecological changes that have consequences for carbon cycling (Chapin, 2005;Hinzman et al., 2005;McGuire et al., 2009;Serreze et al., 2000;Smith and Dukes, 2012;Walther, 2010;Wu et al., 2012).Annual mean surface air temperatures over land increased by Published by Copernicus Publications on behalf of the European Geosciences Union.L. R. Welp et al.: Increasing summer net CO 2 uptake in high northern ecosystems 0.64 • C decade −1 north of 60 • N from 1979 to 2008, roughly twice the rate of 0.33 • C decade −1 for the Northern Hemisphere as a whole (Bekryaev et al., 2010).This northern polar amplification has been attributed to ice/snow-albedo feedbacks (Cess et al., 1991;Qu and Hall, 2007;Serreze and Barry, 2011).Minimum sea ice extent in the Arctic Ocean has declined rapidly (Comiso et al., 2008), with feedbacks and teleconnections on the continental areas as well (2016;Francis et al., 2009).Impacts in the northern regions are predicted to intensify, as climate scenario modeling projects further arctic temperature increases of 5-7 • C by the end of this century (ACIA, 2004;Ewers et al., 2005), and atmospheric CO 2 concentrations continue to rise at a rate approaching 2.0 ppm per year.
Tundra ecosystems and boreal forests hold large stores of carbon in soil organic matter buried in cold or frozen permafrost soils.It is estimated that 1400 to 1850 Pg C are stored in high northern latitude soils and another 60 to 70 Pg C in above-and belowground vegetation (Gedney et al., 2006;Keenan et al., 2013;McGuire et al., 2009).The natural turnover time of this carbon is very slow, but there is a risk that warmer temperatures will increase microbial respiration rates and expose previously frozen organic matter to decomposition by melting the permafrost (Cao et al., 2010;Johnstone et al., 2010;Schuur et al., 2008;Trahan and Schubert, 2015).Over the past several decades there has been a measurable trend to earlier spring snowmelt and surface soil thaw (Brown et al., 2013;McDonald et al., 2004;Smith, 2004) and increases in permafrost borehole temperatures (Romanovsky et al., 2010), demonstrating changes in the thermal stability of northern circumpolar soils.It has been speculated that warming could trigger a massive release of carbon from these soils, in the form of CO 2 and CH 4 , leading to a positive climate-carbon feedback (Pastick et al., 2015;Schuur et al., 2008).This has been nicknamed the arctic "carbon bomb" in the popular media (Treat and Frolking, 2013).
However, warming and the associated lengthening of the growing season encourages plant growth in these otherwise temperature-limited areas, as does increased atmospheric CO 2 fertilization (Lloyd and Farquhar, 1996;Wickland et al., 2006;Yi et al., 2009) and increased nitrogen deposition (Baird et al., 2012;Barber et al., 2000;Bunn and Goetz, 2006;Goetz et al., 2005;Holland et al., 1997;Soja et al., 2007).The net carbon balance of increased plant growth and increased soil respiration is unclear but has important consequences for predicting carbon-climate feedbacks (Abbott et al., 2016;Koven et al., 2015).
Measurements of atmospheric CO 2 concentrations at Barrow, Alaska, by Keeling et al. (1996) provided evidence for increased photosynthetic activity and net primary production (NPP) at northern latitudes from 1960 through 1994.Changes in the phase and amplitude of the CO 2 seasonal cycle were attributed to increased CO 2 uptake by vegetation during spring and summer, leading to earlier drawdown and larger seasonal amplitudes of atmospheric CO 2 concentra-tions (Barichivich et al., 2014;Keeling et al., 1996;Randerson et al., 1999).This perspective is also supported by satellite observations of an increase in vegetation greenness at northern latitudes (Myneni et al., 1997) and global ecosystem process models suggesting that northern ecosystems have become more productive as a result of combined changes in temperature, CO 2 concentration, and nitrogen availability (Kimball et al., 2007;McGuire et al., 2001;Wang et al., 2012).An updated perspective on the northern CO 2 seasonal cycles from Barrow data and from repeated airborne surveys of the mid-troposphere showed a 50 % increase in the amplitude from 1960 to 2010, implying a significant increase in northern ecosystem growing season CO 2 uptake over the last several decades (Graven et al., 2013;Ueyama et al., 2014).
Since the late 1990s, however, some indicators of ecosystem function suggest that the terrestrial biosphere response to recent climate change in the high northern latitudes may be different from the previous few decades and that terrestrial CO 2 uptake has since slowed down or even turned to a net source.Analysis of the changes in the seasonality of atmospheric CO 2 suggests that temperature-induced late summer drought may be increasing fall CO 2 release and offsetting enhanced spring CO 2 uptake (Angert et al., 2005;Piao et al., 2008;Rozendaal et al., 2009).Piao et al. (2008) estimated that current warming during autumn increases respiration in northern ecosystems enough to cancel 90 % of the increased spring CO 2 uptake.While these studies provide important insights into changing ecosystem function, changes in CO 2 seasonal cycles in the atmosphere depend not just on surface fluxes but also on variations in atmospheric circulation (Higuchi, 2002).
Vegetation productivity and distribution have also changed during this same period.The tree line has advanced northward and woody shrub colonies have expanded in the tundra zone displacing less productive species (Goetz et al., 2005;Lloyd et al., 2005;Pearson et al., 2013;Sturm, 2005;Tape et al., 2006).Within the boreal zone, satellite observations show that large areas of the boreal forest not disturbed by fire have been "browning" since 2000 as observed by normalized difference vegetation index (NDVI) measurements (Goetz et al., 2005(Goetz et al., , 2007;;Verbyla, 2008Verbyla, , 2011;;Zhang et al., 2008).These studies are mostly based on either trends in maximum or growing season integrated NDVI or NDVI-derived leaf area index (LAI).One statistical analysis suggests that the browning trends in the Alaskan boreal forest have been ongoing for the last 3 decades (Forkel et al., 2013).These results have been consistent with ground observations which also report widespread tree mortality caused by insect outbreaks due to warmer winter temperatures (Kurz et al., 2008) and drought (Hogg et al., 2008;Peng et al., 2011).However, an updated NDVI processing algorithm in the GIMMS 3g product shows overall more areas greening than the older version, with the largest greening in western Eurasia (Xu et al., 2013;Zhu et al., 2016) but with browning still occurring in parts of North America (Xu et al., 2013).Park et al. (2016) estimates that the overall greening, or increase in growing season integrated NDVI, is equivalent to a 20.9 % gain in productivity since 1982 and the smaller areas of browning are equivalent to a 1.2 % loss of productivity.Furthermore, these trends in greening and productivity seem to be independent of shifts in the start and end of the growing season and growing season length (Park et al., 2016).
It is important to determine the net carbon balance of these large northern regions to see if they have been increasing or decreasing in CO 2 sink strength or perhaps transitioning to a net CO 2 source.Approaches used to estimate the net CO 2 fluxes of large areas include forest inventories, atmospheric inversions, and process-based models.Each of these methods has its strengths and weaknesses.Atmospheric inversions can infer the global, continental, and sometimes regional-scale fluxes of CO 2 between the atmosphere and the land biosphere and the oceans by analyzing the temporal and spatial records of atmospheric CO 2 change (Enting, 2002).Inversions have the advantage of including the effects of disturbance and changing vegetation patterns but are limited to the period of sufficient CO 2 concentration observations and are best suited to resolving continental-scale fluxes.Few inversion analyses have specifically focused on the high northern latitude terrestrial ecosystems.Zhang et al. (2013) aggregated the inversion fluxes from five different models into Eastern Canada and Western Canada plus Alaska.The inversions examined showed consistent increases in CO 2 uptake in Eastern Canada and no long-term trend in Western Canada plus Alaska.
Atmospheric inversion CO 2 fluxes can provide useful validation metrics for other methods of ecosystem monitoring because they resolve the net effect of above-and belowground carbon fluxes.Repeat forest inventories are useful for identifying trends in forest productivity over time.However, they are limited to detecting trends in aboveground carbon only, not changes in belowground carbon.Forest inventory studies have found drought-induced tree mortality and aboveground carbon loss in Canada, with the western region most affected (Hogg et al., 2008;Ma et al., 2012;Michaelian et al., 2010;Peng et al., 2011).
Process-based modeling studies attempt to account for above-and belowground carbon changes while providing full spatial coverage and are therefore capable of simulating net ecosystem fluxes that can be compared to atmospheric inversions.Several process-based modeling studies have concluded that the arctic tundra and boreal forests have been decreasing sinks or increasing sources since the 1980s due to climate effects, namely warmer temperatures increasing soil organic matter decomposition, and increased fire and insect disturbance, offsetting increased CO 2 uptake driven by CO 2 fertilization (Bradshaw and Warkentin, 2015;Hayes et al., 2011;McGuire et al., 2010), though McGuire et al. (2012) found that 6 out of 11 process-based models (global and regional) predicted a strengthening carbon sink (including CH 4 contributions) in the 2000s compared to the 1990s.Pro-cesses including permafrost melt, hydrologic changes, nutrient dynamics, and fire emissions are critical to predicting any changes in the net CO 2 fluxes from the northern regions but difficult to include in the models with much certainty (Abbott et al., 2016;Harden et al., 2012).Abbott et al. (2016) conducted a survey of experts in the field and found that the overwhelming opinion was that any increases in biomass are not going to be enough to offset carbon releases from permafrost melt by 2040 or 2100.These models and expert predictions can be directly challenged with results from atmospheric inversions.
In this study, we examine large regional-scale temporal variability in terrestrial CO 2 fluxes from two atmospheric inversions using interannually variable atmospheric transport from 1985 to 2012.We focus on trends in the carbon uptake of the land biosphere north of approximately 50 • N. The primary objective of this study is to evaluate temporal changes in the annual and seasonal land biosphere CO 2 fluxes.We determine in what months surface CO 2 fluxes have likely changed, i.e., increased summer uptake or winter release.We further examine NDVI, air temperature trends, and correlations with CO 2 fluxes to provide some spatial and process context for temporal changes in the inversion fluxes.
2 Methods and data analysis

Atmospheric inversions
We compared two different atmospheric inversions: the RIGC inversion and the Jena CO 2 inversion (s85 v3.6).The RIGC inversion method was adapted from Rayner et al. (1999), and largely followed the TransCom-3 protocol (Gurney et al., 2003).The RIGC inversion uses a 64-region time-dependent inverse method to infer carbon source and sink estimates based on the method of Patra et al. (2005).The RIGC inverse calculation starts with a priori fossil fuel emissions and terrestrial and oceanic fluxes, which are then optimized to match observations.For the a priori fluxes, terrestrial exchanges were taken from the Carnegie-Ames-Stanford Approach (CASA) model monthly output (Randerson et al., 1997) and monthly mean oceanic exchanges from Takahashi et al. (2002) as in the TransCom-3 protocol (Gurney et al., 2003).Total anthropogenic CO 2 emissions were derived from the Oak Ridge National Lab monthly fossil fuel estimates from the Carbon Dioxide Information Analysis Center (CDIAC) (Boden et al., 2009) plus bunker fuel and non-fuel oxidation estimates from the Emissions Database for Global Atmospheric Research (EDGAR) (Oliver and Berdowski, 2001) derive residual a posteriori land and ocean surface fluxes for the 64 inversion regions.
In the RIGC inversion, the GLOBALVIEW-CO 2 input was limited to the 26 stations, which have nearly continuous CO 2 observation from 1985 to 2006 (Table 1).Among these are stations that document changes in high northern latitudes, including Barrow (71 • N), Alert (82 • N), Station M (66 • N), Cold Bay Alaska (55 • N), Shemya (52 • N), and Cimone (44 • N).A selected set of stations was used to avoid creating spurious trends in the inversion results from adding new stations midway through the inversion period.All selected stations had at least 71 % of months sampled at least once and came online by 1989.Stations north of 39 • N had 84 to 100 % monthly coverage.The resulting fluxes from the RIGC inversion are valid from 1986 to 2006, after removing years at the beginning and end for "edge effects" from the inversion setup.
Compared to the RIGC inversion, the Jena inversion, run s85 v3.6 (Rodenbeck, 2005), uses a slightly different set of 19 stations selected to completely cover the 1985-2012 estimation period but includes all of the same stations north of 50 • N (Table 1).It uses individual measurements from various sampling networks, without smoothing or gap filling.Fluxes are estimated on a grid-scale resolution (approximately 4 • latitude by 5 • longitude) to reduce aggregation errors.However, to avoid underdetermined solutions, spatial and temporal a priori correlations are imposed, smoothing the estimated flux field on scales smaller than about 1 week and about 1600 km (land, in longitude direction), 800 km (land, latitude), 1900 km (ocean, longitude), or 950 km (ocean, latitude).Land flux adjustments are spatially weighted with a productivity proxy, the long-term mean NPP from the Lund-Potsdam-Jena (LPJ) terrestrial biospheric model (Sitch et al., 2003).Prior fluxes comprise anthropogenic CO 2 emissions, including fossil and bunker fuels, from EDGAR v4.2 (EDGAR, 2011), a constant spatial flux pattern on land (time mean net ecosystem exchange (NEE) from the LPJ model), and an ocean-interior inversion by Mikaloff Fletcher et al. (2006), with a mean seasonal cycle of ocean fluxes from Takahashi et al. (2002).The Jena inversion uses the TM3 global atmospheric transport model driven by interannually varying meteorology from the NCEP reanalysis.The 4 × 5 • gridded a posteriori land and ocean surface fluxes are aggregated to our larger analysis regions.Resulting fluxes are valid from 1985 to 2012.

Datasets
We compared CO 2 fluxes with satellite-based NDVI data over the same time period and with land temperature records.NDVI is a proxy for photosynthetically active aboveground biomass, calculated from the visible and near-infrared light reflected by vegetation.It has dimensionless units and varies from a value of 0 for no vegetation to a value of 1 for the highest density of green leaves.We used NDVI data produced by NASA's Global Inventory Modeling and Mapping Studies (GIMMS 3g) from measurements of the Advanced Very High Resolution Radiometer (AVHRR) satellite and supplied at the monthly, 1 × 1 • resolution (Pinzon and Tucker, 2014).Winter NDVI data were excluded from this analysis because of the confounding influence of snow (Myneni et al., 1997).We defined growing season NDVI (Zhou et al., 2001) as the sum of monthly NDVI from April to October following the example of earlier work.Comparisons with recent satellite measurements of solar-induced fluorescence show that the seasonality of NDVI may not fully capture the seasonality in gross primary production (GPP) (Walther et al., 2015), but we focus on interannual variability of growing season sums and maximum July values in this study.
We examined estimates of fire CO 2 emissions from 1985 to 2000 from the RETRO reanalysis project and from 1997 to 2012 from the Global Fire Emissions Database (GFED4) (Giglio et al., 2013;Schultz et al., 2008).

Analysis approach
Our focus is primarily on the interannual variability in the CO 2 fluxes, which is considered more robust across inversions than the absolute values of the mean fluxes (Baker et al., 2006).For that reason, we look at anomalies from the long-term mean values of each inversion.
We focused our analysis on land carbon fluxes in two roughly zonal bands at high northern latitudes partly based on the TransCom regional boundaries defined by Gurney et al. (2003).Figure 1 shows the regions of boreal Asia (BA) and boreal North America (BNA) that we aggregated into what we refer to as the "boreal zone", roughly between 50 and 60 • N, and the "arctic zone", north of 60 • N. Note here that while we refer to the boreal zone as roughly 50 to 60 • N, the southern boundary is not defined at the 50 • N latitude but follows the irregular southern boundary of boreal forest (stippled area in Fig. 1) nor does it include the entire boreal forest as the northern boundary extends well into the arctic zone north of 60 • N. We decided to omit the European (EU) land region from our zonal analysis for two reasons.First, the TransCom protocol followed by the RIGC inversion does not separate northern Europe at 60 • N like it does for BA and BNA; rather, the northern EU section is everything north of 50 • N. Second, the EU region includes a relatively small fraction of the tundra and boreal forest ecosystems compared to BA and BNA, and the forest area is highly managed.Our focus is how the less intensively managed ecosystems of the north have been responding to climate change, and examining the boreal zone and the arctic zone should maximize any potential signals of change.A similar approach of excluding EU was used in the arctic analysis of McGuire et al. (2009).
The atmospheric inversion approach taken in this study is unlikely to reliably separate influences from different longitudinal regions within the latitude bands discussed here.Our focus on the longest records possible, from sparse atmospheric CO 2 observations starting in the 1980s, compromises the spatial resolution of the inversion fluxes.Rapid atmospheric mixing of a few weeks around latitude bands makes it hard to separate fluxes for example from North America and Eurasia.Spatial errors in the assignment of surface fluxes by the inversion analyses could create spurious trends and complicate our interpretation of the data.For that reason, we check the EU and northern ocean (NO) regions for any trends that might be offsetting trends in what we define as the boreal zone and the arctic zone.We performed the trend analysis for two periods: from 1986 to 2006, when we have results from both inversions, and from 1985 to 2012 for just the Jena inversion.
We examine trends in the monthly, seasonal, and annual NEE fluxes from the inversion models.Amplitudes of the annual seasonal cycle in CO 2 fluxes were calculated from the maximum and minimum monthly mean fluxes within each calendar year as flux amplitude = maximum NEE-minimum NEE.The flux amplitude is indirectly related to the ampli- tude in atmospheric CO 2 concentrations, as the atmospheric concentration is roughly the integral of the monthly fluxes.It is unnecessary to detrend the time series of fluxes from the inversions prior to calculating the flux amplitude, unlike the concentration amplitude, which has a persistent long-term trend from anthropogenic CO 2 emissions.We also examine the latitudinal gradient of the trends in the seasonal fluxes in ∼ 4 • latitude bands from the gridded Jena inversion.This analysis was not possible with the RIGC inversion because of the larger basis regions.This approach attempted to answer the question of whether summer uptake or fall respiration (or both) is increasing and how that might change with latitude.Two statistical methods were used to calculate the slopes and significance of long-term trends: linear least squares (LSQ) and Mann-Kendall (M-K) tests.Trends were considered significant if they passed the 90 % confidence level  (p values < 0.1).The Model I linear region analysis (LSQ) was done in Matlab using the function "lsqfity" developed by Peltzer (2000) based on Bevington and Robinson (1992).
The nonparametric monotonic Mann-Kendall trend test with Sen's slope was also done in Matlab using the function "ktaub" developed by Burkey (2006).Results of both tests are presented in Table 2.The independent tests generally agree on slope and significance of trends.2).Negative values represent uptake of CO 2 by the land biosphere, i.e., out of the atmosphere.The arctic zone containing the tundra region showed no significant trend in annual CO 2 uptake (Fig. 2a, Table 2) from 1986 to 2006 in either inversion.The longer period, 1985-2012, in the Jena inversion did show a small but significant trend toward increased uptake of an extra 4 Tg C yr −1 .Anomalously strong annual CO 2 uptake occurred in the years 1990 and 2004 in the RIGC inversion and strong CO 2 release in 1996.These large anomalous fluxes were not present in the Jena inversion.The Jena and RIGC inversions differ in their mean seasonal cycle in the arctic zone, with the RIGC inversion yielding peak CO 2 uptake approximately twice that of the Jena inversion (Fig. 3b).The seasonal amplitude and phase has previously been found to differ among atmospheric inversions (McGuire et al., 2012).These differences are not unexpected given the differences in atmospheric transport (including vertical mixing and leakage across latitudes), a priori fluxes, observational network inputs, and model structure between the inversion.In this analysis we try to focus on the most robust features where the inversions do tend to agree on the interannual trends in anomalies from the mean.Trends in monthly net CO 2 flux, computed with the method of Randerson et al. (1997), reveal increasing uptake in July in both inversions and stronger releases in September, October, and November (Fig. 3d).These seasonal changes largely cancel in the annual net fluxes but contribute to increasing CO 2 flux amplitudes, computed as the difference between the maximum and minimum monthly CO 2 fluxes, by ∼ 1.0 % yr −1 relative to the mean seasonal amplitude from 1986 to 2006 for both inversions (Fig. 4, Table 2).Figure 5 shows the annual values of the July CO 2 flux in Pg C yr −1 over this record.This is directly related to July trend data in Fig. 3d.On a per area basis, this translates into an increase in July peak summer CO 2 uptake of 0.007-0.013gC m −2 d −1 yr −1 , depending on the inversion used, averaged over the entire zone or a ∼ 10 % increase in peak summer CO 2 uptake over these 21 years.The trends over 1985-2012 are similar at 0.007 gC m −2 d −1 yr −1 for the Jena inversion (Table 2).

Boreal zone (50-60 • N)
The boreal zone shows a trend towards increasing annual net CO 2 uptake in both inversions (Fig. 2b, Table 2).From 1986 to 2006, the trend in the RIGC inversion was 10 Tg C yr −2 with a p value < 0.1.The Jena inversion resulted in a similar trend of 8 Tg C yr −2 but did not meet the criteria for significance, p > 0.1 (Table 2).The most noticeable difference between the inversions is that the RIGC inversion predicted an anomalous release of CO 2 in 1994 that was not confirmed by the Jena inversion.Over the longer period from 1985 to 2012, the Jena inversion predicts the same trend toward greater CO 2 uptake, with a slope of 7-8 Tg C yr −2 and a p value < 0.1 (Fig. 2b, Table 3).
The Jena and RIGC inversions resulted in similar mean seasonal cycles of the monthly net CO 2 fluxes in the boreal zone, but the seasonal amplitude in the Jena inversion was slightly larger (Fig. 3a).The trends in the monthly fluxes show increasing CO 2 uptake in the growing season, and in the case of the Jena inversion, increasing CO 2 uptake in the spring and release in the fall (Fig. 3c).There was a corresponding increase in the seasonal amplitude of net CO 2 flux of 0.4 % yr −1 (p = 0.04) estimated by the Jena inversion but not the RIGC inversion (Fig. 4b and Table 2).Figure 5 shows the time series of CO 2 flux in July (month of peak flux) over this period.Both inversions show an increase in July CO 2 uptake although they do not agree on anomalies from year to year.In Fig. 3c, d, both inversions also show an increase in the fall CO 2 release in the northern land regions, but the Jena inversion attributes this mostly to the boreal zone, whereas the RIGC inversion attributes it mostly to the arctic zone.

Europe and northern ocean fluxes and fossil fuel emissions
For completeness, we also show the time series of CO 2 flux trends, both net annual and seasonal amplitude from the 55 to 80 • N region of EU and the NO, to be sure that fluxes in these regions are not compensating for fluxes in our analysis of the BA + BNA regions (Fig. 6).There were no offsetting positive trends in the annual net flux of CO 2 or negative trends in the seasonal amplitude in EU from 1986-2006 and none were statistically significant (Table 2).Likewise the NO flux trends are insignificant with the exception of the seasonal amplitude trend in the RIGC inversion of −2.4 % yr −1 .This was statistically significant using the modified Mann-Kendall p test, but the small percentage increase in the mean amplitude of the ocean flux (∼ 0.1 Pg C yr −1 ) is much too small to offset gains in the land flux amplitude in the BA + BNA regions (mean 4.6-9.9Pg C yr −1 for north of 60 • N and mean 11.4-13.9Pg C yr −1 for 50-60 • N).There is no indication that inversion-resolved trends in the EU and NO regions in the north of the 50 • N zone are forcing offsetting trends in the BA + BNA regions; however, we cannot   in absolute fluxes were similar in each inversion.Differences in the fossil emissions are therefore unlikely to contribute significantly to trends in the biological land fluxes of the BA + BNA arctic and boreal zones.

Flux amplitude trends
We define the seasonal flux amplitude as the difference between the peak summer CO 2 uptake and the maximum CO 2 release in the fall, within a calendar year.We examined changes in the flux amplitude using several approaches.Figure 4 shows that the flux amplitude increase, in percent of the mean flux amplitude, is larger in the arctic zone than the boreal zone.This is also reflected in the monthly trends in Fig. 3c, d.
We also examined the change in the seasonal flux amplitude across latitudes from the Jena inversion to see if this observed increase in the seasonal flux amplitude was unique to the high northern latitudes or if it is more widespread.Here we define the fall flux as the mean of SON, and the summer uptake is fixed as July.Figure 7a shows that significant increases in the annual flux amplitude have occurred between 40 and 70 • N with a peak from 50 to 65 • N. Looking at the fall and summer contributions separately shows that both increasing fall CO 2 release and peak summer uptake contribute to the annual amplitude increase, with increasing fall CO 2 release outpacing peak summer uptake in the 40-50 • N band and increasing peak summer CO 2 uptake outpacing fall release in the 55-70 • N band (Fig. 7b).

Fire emissions
The net fluxes examined here are dominated by land biosphere fluxes, but they also include CO 2 emissions from forest fires, which occur mostly in summer (Fig. S3 in the Supplement) (van der Werf et al., 2006).If fire activity were responsible for the trend in summer net carbon uptake, then fire emissions would need to be decreasing.We examined estimates of fire CO 2 emissions from 1985 to 2000 from the RETRO compilation and from 1997 to 2012 from the GFEDv4 model (Giglio et al., 2013;Schultz et al., 2008).While we cannot combine the two emissions estimates into a continuous time series because of the different methodologies used, we can examine trends over each record.Neither record shows evidence of decreasing fire emissions over their respective time periods (Fig. S4 in the Supplement); therefore, biological activity clearly dominates the trends.

Temperature and NDVI trends
In order to investigate possible drivers of the trends in CO 2 fluxes, we also examine trends in surface air temperature and NDVI.Seasonal temperature changes from 1986 through 2006 were not uniform across the far north (Fig. 8).In general, warming has been the greatest in the fall (SON) and winter (DJF), although these patterns vary regionally.Despite the general trend toward warming, cooling trends are seen over the 1986-2006 period for Siberia in winter (DJF) and western North America in spring (MAM).Nearly all of the land regions in the Northern Hemisphere have experienced warmer summers (JJA) and falls (SON).
Figure 9 shows linear trends from 1986 to 2006 in gridded NDVI averaged over the "growing season" (April through October) and for the month of July.Widespread greening trends are observed, with the exception of browning in the southern boreal forest of North America.Significant greening trends are found in tundra regions, especially in North America.
Figure 10a averages the growing season NDVI over the latitude bands, for each year, showing both the growing season average (top panel) and seasonal maximum (bottom panel).The period 1986-2006 showed no significant trend in the 50-60 • N band in either growing season or maximum metrics (Fig. 9).In contrast, a significant increasing trend of 0.13 % yr −1 (Table 2, p = 0.01) is found in the > 60 • N band over this period, driven mostly by trends in tundra regions.Although not included in the trend analysis, growing season NDVI north of 60 • N increased abruptly at the end of the record by ∼ 5 % from 2009 to 2010.Similarly large changes occurred in 1991-1992 and 1996-1997.Figure 10b shows the trend in annual peak NDVI, the maximum monthly value for each year regardless of which month it is.This also shows a small but significant increase of 0.12 % yr −1 (Table 2, p = 0.0103) north of 60 • N from 1986 to 2006 and no significant trend in the 50-60 • N region.Compared to growing season NDVI, the peak NDVI increase from 2009 to 2010 was less extreme.Plant growth in 2010 was increased in the shoulder seasons in addition to the mid-season peak.

Controls of temperature and NDVI on CO 2 fluxes
Summer uptake and fall release of CO 2 play a large role in the atmospheric CO 2 concentration, so here we look at the correlations of CO 2 fluxes with temperature and NDVI as  proxies for soil respiration and primary productivity variability to help assess mechanistic links.We performed lagged correlation analysis on monthly time series by calculating the temporal correlation for either the July net CO 2 fluxes or fall   (SON) fluxes and 3-month running means of temperature or NDVI time series with 0-to 60-month (up to 5-year) lags.In this analysis, all data sets were de-trended using a stiff spline to remove long-term trends, thus emphasizing processes controlling interannual variability.It is also important to remember that NEE is negative when there is net CO 2 uptake from the atmosphere when interpreting the sign of correlations.
Figure 11 shows the lagged correlations for the period 1986-2006 of both inversions for the boreal zone.We found significantly strong positive correlations between July CO 2 NEE and April through August temperatures of the same year but no evidence of correlation with NDVI.The temperature correlation suggests that warmer growing season temperatures increase soil respiration (or wildfires) and result in reduced peak CO 2 uptake (less negative NEE).The RIGC inversion shows a significant correlation between warm winters and increased CO 2 uptake the following growing season (negative correlation), consistent with Russell and Wallace (2004), but this relation did not appear in the Jena correlation.The fall CO 2 fluxes were also positively correlated with growing season temperatures of the same year (not shown).This is consistent with increased growing season air temperatures stimulating soil respiration through the fall, either through increased carbon pools from enhanced summer productivity or warmer soil temperatures that persist into the fall.
Our analysis found significant correlations up to 2 to 4 years before, suggesting that temperature anomalies could have an impact on NEE after several years' delay.While there have been studies suggesting that multiyear lags between climate and CO 2 fluxes are important (e.g., Bond-Lamberty et al., 2012), these correlations were not consistent between the inversions, preventing us from speculating as to the cause.
An analysis of the peak July NDVI correlations with lagged air temperatures showed that warmer temperatures in the 14 months prior to the peak NDVI were associated with higher NDVI in the boreal zone (Fig. 12).One possible interpretation of the correlation analyses presented here is that warmer temperatures in the boreal zone lead to increased plant productivity (indicated by positive NDVI and temperature correlation) but that the IAV of the net C balance in July and the fall was dominated by respiration (indicated by positive NEE and temperature correlation and by the lack of an NDVI and NEE correlation).
Lagged correlations for the arctic zone were generally less significant (Fig. S2 in the Supplement).July CO 2 fluxes did not show a strong correlation with either temperature or NDVI.Fall CO 2 fluxes were not consistently correlated with temperatures in the same season in either inversion.Peak July NDVI was weakly correlated with May-June temperatures.Overall, we conclude that this northern region may be dominated by other controls, like soil thermal processes that would not show up clearly in the correlation analysis with air temperature and NDVI.

Discussion
The results of these two inversions show that the northern high-latitude regions of BNA+BA remain nearly constant or slightly increasing sinks of atmospheric CO 2 .The boreal zone, again excluding Europe, absorbed an extra 8-11 Tg C yr −2 over the period from 1986 to 2006, resulting in an annual CO 2 sink in 2006 that was 170-230 Tg C larger than in 1986.This trend towards increasing CO 2 uptake appears to continue through 2012 as indicated by the longer Jena inversion.This result contradicts some modeling stud-ies, which point to trend reversals in observed NDVI and modeled net CO 2 fluxes.Hayes et al. (2011) used the TEM ecosystem model to show that increased respiration and fires in this region had weakened the sink strength since 1997.Dynamic global vegetation models have also predicted a trend toward CO 2 release to the atmosphere across much of the northern land region from 1990 to 2009 (Sitch et al., 2015).In general, Carvalhais et al. (2014) found that models tend to overpredict the transfer of carbon from the soils to the atmosphere and overestimate the sensitivity of heterotrophic respiration to climate, although they did not include the TEM model used by Hayes et al. (2011).The results of this inversion analysis suggest that respiration is overpredicted in models or that increased primary production not captured by these models is offsetting increases in soil respiration and/or forest fire emissions.Our results are also consistent with Forkel et al. (2016), who use process modeling constrained by atmospheric observations to conclude that photosynthesis has responded more strongly to warming than carbon release processes.
Sensor drift and calibration errors may have resulted in false browning trends in the boreal forest, particularly in the needleleaf evergreen forests.Recent analyses of NDVI trends in the updated GIMMS 3g version find significantly more greening trends than in the previous GIMMS version (Bi et al., 2013;Guay et al., 2014).On a pan-boreal basis, it seems plausible that the CO 2 sink strength has continued to increase despite previous reports of drought stress reducing CO 2 uptake of the boreal region.As a complication, however, any changes in net carbon fluxes in these ecosystems will depend not only on aboveground vegetation changes that can be observed by remote sensing but also on processes occurring below ground, where most of the carbon is stored (Iversen et al., 2015).
For the arctic zone, we estimate that July CO 2 uptake increased from 1986 to 2006 by 0.15 to 0.27 g C m −2 d −1 , depending on the inversion and the trend detection algorithm.This estimate is based on multiplying the regression slopes (Table 2) by the 21-year time frame.In this zone, we found the strongest NDVI greening trends in the tundra regions, covering roughly 25 % of the relevant land area.In light of evidence of rapid shrub expansion in these tundra ecosystems (Myers-Smith et al., 2011;Tape et al., 2006) (Elmendorf et al., 2012), a rapid increase in July CO 2 uptake by the tundra ecosystems is plausible.
Most of the previous studies investigating seasonal variability in the northern ecosystem carbon fluxes have relied on observations of atmospheric CO 2 concentrations (e.g., Angert et al., 2005;Buermann et al., 2013;Keeling et al., 1996;Piao et al., 2008;Randerson et al., 1999).The analysis presented here is unique in that it considers variability in atmospheric transport through the atmospheric inversion approach.Our results are generally consistent with the finding of Piao et al. (2008) that enhanced CO 2 losses from northern ecosystems in the fall partially cancel the enhanced CO 2 uptake earlier in the growing season, especially in the arctic zone.We also find evidence of uptake enhancement in the summer as well as the spring in both the arctic and boreal zones, consistent with Graven et al. (2013) and Forkel et al. (2016).
Our investigation of the controls on interannual variability of the CO 2 fluxes showed increased CO 2 uptake in cooler summers.There are many previous attempts to identify the short-term drivers of the net carbon balance of ecosystems from eddy covariance studies.Several studies have found that warm and dry summers lead to drought stress and reduced net CO 2 uptake in boreal forest ecosystems (Arain et al., 2002;McMillan et al., 2008;Welp et al., 2007).Net CO 2 flux reductions could be the result of decreased primary productivity, increased respiration, or both.Correlations between temperature and annual tree ring growth increments point to a switch from a positive correlation to a negative correlation (reduced growth during warm years) driven by increased drought stress in recent decades (Barber et al., 2000;Beck et al., 2011). Wunch et al. (2013) and Schneising et al. (2014) found that summertime total-column CO 2 in the north was relatively higher during years with warm anomalies in the boreal region, suggesting that reduced net CO 2 uptake during warm summers may be driven by the temperature dependence of soil respiration.These correlations of interannual variability are contrary to the overall long-term association with warming and greater summer uptake (Keeling et al., 1996).This difference could reflect the importance of structural ecosystem changes due to warming on a long timescale, increasing photosynthesis (Graven et al., 2013), but are also consistent with respiration as the dominant control of NEE on short timescales (Schaefer et al., 2002).The difference in the timescale of response is important to consider.
A full explanation of the trends in CO 2 fluxes of the arctic and boreal zones is still lacking, with possible causes including changes in temperature, which were explored here but also soil moisture, nutrient status, or fire and insect disturbance.An important unresolved question is how the distribution of deciduous and evergreen plant functional types has changed on the pan-boreal scale over this period.A shift to younger forests, with an increasing deciduous fraction, would increase the seasonal flux amplitude (Welp et al., 2006;Zimov, 1999), perhaps with little change in common NDVI metrics.
The latitude gradient of changes in land fluxes from the Jena inversion, now including results from Europe, shows that, from the 1985-1989 mean to the 2007-2011 mean, the surface flux amplitudes have increased the most in the zone from 40 to 65 • N (Fig. 7a).This is consistent with Graven et al. (2013), who argued that, over the longer period from 1960 to 2010, the increases in CO 2 flux amplitude were centered mostly on boreal regions.Our analysis of the Jena inversion by latitudinal bands shows that the increases in peak July uptake have been greater than the fall CO 2 releases north of 55 • N, but from 40 to 50 • N, fall release outpaced July up-take (Fig. 7b).The results presented here show that the increased seasonal amplitude in atmospheric CO 2 in the high northern latitudes is not caused by flux trends in the summer or fall only but rather that both contribute (Fung, 2013;Graven et al., 2013).The trend in annual net CO 2 fluxes also includes changes in other months, namely greater uptake in the spring (not shown in Fig. 7), which contributes to the annual sum.The advantage of this analysis is that it incorporates interannually varying atmospheric transport, so the temporal changes in the surface fluxes should be better resolved.It does not identify whether individual months have a disproportionately larger influence on the atmospheric CO 2 concentration amplitude.
Our attempt to distinguish changes by ∼ 10 • latitude arctic and boreal zones is pushing the limits of what is feasible from atmospheric inversions based on sparse atmospheric CO 2 observations.The limitation is illustrated by the tendency of the two inverse calculations to allocate the increase in fall CO 2 release mostly to different latitude bands and the shift to increasing earlier CO 2 uptake in the Jena inversion compared to the RIGC inversion in Fig. 3c.Resolving fluxes with a monthly resolution is also challenging (Broquet et al., 2013), but the long record examined here, by two independent inversions, gives us reasonable confidence in this aspect.

Conclusions
The two atmospheric inversions analyzed in this study show that the annual net CO 2 sink strength in the boreal zone has increased from 1985 to 2012.However, the annual net CO 2 fluxes in the arctic zone showed no trend.Both regions show significant increases in midsummer CO 2 uptake.But a trend towards greater CO 2 emissions in the fall has partly canceled the trend toward greater summer uptake, with the largest cancelation in the arctic zone.These trends in summer and fall fluxes cause the seasonal amplitude of the fluxes to increase and, consequently, the seasonal amplitude of atmospheric CO 2 concentrations.
We also examined NDVI and CO 2 flux interannual correlations showing that while warmer summers were correlated with increasing NDVI in the long term, relatively cooler summers favor net CO 2 uptake in the boreal region in the short term.This suggests that increased respiration can outpace increases in productivity in the short term.Overall, there is evidence from these atmospheric inversions that increased CO 2 uptake from the northern region is offsetting carbon release in the pockets of browning in the boreal zone.Our findings are consistent with the recent NDVI studies from the GIMMS 3g product that find overall greening of the boreal and arctic regions (Park et al., 2016;Xu et al., 2013;Zhu et al., 2016).By itself, this would suggest the potential for continued or strengthening net CO 2 uptake.However, north of 60 • N, our findings show that fall CO 2 release largely offsets increased summer uptake in the net annual budget.These results un-derscore the difficulty of resolving net fluxes from remotesensing indices alone, which can only "see" the productivity and not the respiration fluxes.
Furthermore, our atmospheric inversions results show no evidence of an overall trend towards increasing CO 2 releases in either the boreal or arctic zone over the 1985-2012 period.This is an important check for process-based biospheric models, which tend to predict a shift from sink to source in the next century (Treat and Frolking, 2013).Our results are not consistent with studies that suggest carbon sinks have weakened in boreal and arctic ecosystems over past decades (Bradshaw and Warkentin, 2015;Hayes et al., 2011) but support some process models which do predict a strengthening CO 2 sink in the northern region (McGuire et al., 2012).To date, the increase in biomass productivity has appeared to be outpacing CO 2 losses from warming northern carbon-rich soils.Time will tell whether this trend continues or whether it will reverse and a net carbon source will develop in a few decades as predicted by popular opinion among the community of experts (Abbott et al., 2016).

Data availability
The data used in this analysis are publically available from the individual authors responsible for creating the products.The Jena CO 2 inversion results are posted on the project website: http://www.bgc-jena.mpg.de/~christian.roedenbeck/download-CO2/.Run ID s85 version 3.6 was used in this project.Associated files contain the atmospheric monitoring site locations and data used in the inversion and the fossil fuel emissions that were used to solve for the biological land CO 2 fluxes.The RIGC CO 2 inversion results are posted on the Global Carbon Atlas project website: http://www.globalcarbonatlas.org/?q=en/ content/atmospheric-inversions.Likewise, associated files contain the atmospheric monitoring site locations and data used in the inversion and the fossil fuel emissions that were used to solve for the biological land CO 2 fluxes.The GIMMS 3g NDVI data used is posted on the AVHRR website, https://nex.nasa.gov/nex/projects/1349/.The Matlab and GrADS scripts used in this analysis are publically available on GitHub: https://github.com/lwelp/Welp_2016_ACP.
The Supplement related to this article is available online at doi:10.5194/acp-16-9047-2016-supplement.
the NSF under award PLR-1304270, and the UC Multiple Campus Award Number UCSCMCA-14-015.Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NASA, the NSF, the DOE, or UC.Prabir K. Patra was partially supported by MEXT Arctic GRENE (ID 5).
Edited by: P. O. Wennberg Reviewed by: three anonymous referees

Figure 1 .
Figure 1.The major land and ocean basis regions used in the RIGC inversion based on the TransCom3 regions.The Jena inversion was done on a ∼ 4 × 5 • grid and aggregated to these regions.The northernmost land regions are shown in color.The two zones that we discuss in this analysis cover boreal North America and boreal Asia and are marked in shades of blue, with the arctic zone (> 60 • N) in light blue and the boreal zone (50-60 • N) in dark blue.The European basis region, in red, is not divided at 60 • N in the RIGC inversion and therefore is not included in this analysis.Stippling indicates the boreal forest biome based on the Global Land Data Assimilation System (GLADS) University of Maryland modified International Geosphere-Biosphere Programme (IGBP) land classification scheme (http://ldas.gsfc.nasa.gov/gldas/GLDASvegetation.php).The tundra biome is north of the stippling.

Figure 2 .
Figure 2. Annual CO 2 fluxes normalized by subtracting the 1986-2006 mean value for (a) arctic zone (> 60 • N) and (b) boreal zone (50-60 • N).Black shows the RIGC inversion results.Grey shows the Jena s85 inversion results.Dashed lines are linear trends from 1986 to 2006 (value in Table2).Negative values represent uptake of CO 2 by the land biosphere, i.e., out of the atmosphere.

Figure 3 .
Figure 3. Mean monthly CO 2 fluxes for (a) arctic zone (>60 • N) and (b) boreal zone (50-60 • N).Black circles are the 1986-2006 means of the RIGC inversion.Grey squares are the Jena inversion over that same period.Magenta is the Jena inversion average over a longer time period (1985-2012).Differences are likely due to differences in atmospheric transport, including vertical mixing, between the inversions.Linear monthly trends of (c) arctic zone and (d) boreal zone for the same inversions and time periods as in (a) and (b).

Figure 4 .
Figure 4. CO 2 flux amplitude for each year calculated as the maximum monthly flux (positive: CO 2 release to the atmosphere) minus the minimum monthly flux (negative: CO 2 uptake by the biosphere) for (a) arctic zone and (b) boreal zone.Black shows the RIGC inversion; grey shows the Jena s85 inversion.Dashed lines show the linear trends from 1986 to 2006, the common period between the inversions (Table2).

Figure 5 .
Figure 5. July CO 2 flux anomalies for each inversion for the (a) arctic zone and (b) boreal zone normalized by subtracting the 1986-2006 mean value.This is the month of maximum CO 2 uptake in each case (see Fig. 3).The dashed lines are the linear trends from 1986 to 2006, also plotted in Fig. 3c and d.

Figure 6 .
Figure 6.Fluxes from the northern ocean and European basis regions.(a) Annual fluxes and (b) annual flux amplitude for the northern ocean.(c) Annual sum and (d) annual flux amplitude for the EU region.The trends in these fluxes are small and, in the case of EU, in the same direction, compared to the trends resolved for the BA + BNA regions.

Figure 7 .
Figure 7. Latitudinal gradients in the land CO 2 fluxes from the Jena s85 inversion in 5 degree bands.(a) An estimate of the change in the annual CO 2 seasonal flux amplitude (Pg C yr −1 ) by subtracting the July from the SON fluxes at 5 degree latitude bands.Positive values reflect an increase in the peak-to-trough flux amplitude.(b) Green is the difference between the 2007-2011 mean and the 1983-1989mean in the July CO 2 uptake with the sign reversed (here positive is uptake by the biosphere) and magenta is the difference in the mean of September-November fall CO 2 release with conventional sign (positive is release of CO 2 to the atmosphere).

Figure 8 .
Figure8.Gridded temporal trends in surface air temperature from the GISS temperature record (http://data.giss.nasa.gov).Plots were made using software available on the data archive website.

Figure 9 .
Figure 9. Gridded temporal trends in GIMMS 3g NDVI (a) growing season (April-October) mean and (b) July only from 1986-2006.Trends are expressed as percent changes from the mean.

Figure 10 .
Figure 10.Time series of NDVI trends averaged for the analysis regions in this study.(a) Growing season (April-October) mean and (b) annual maximum, usually in July.Black is for the arctic zone and grey is for the boreal zone.

Figure 11 .
Figure 11.Correlation coefficients for July CO 2 fluxes in a given year (year 0) from the boreal zone with lagged 3-month running mean temperature (area-weighted and NPP-weighted) and NDVI for the same region.(a) RIGC inversion and (b) Jena s85 inversion over the current and previous 4 years.Positive correlations mean that high temperature or NDVI leads to less CO 2 uptake.Filled circles indicate significance greater than the 95 % level.Shaded bars indicate the summer months (May-August).

Figure 12 .
Figure 12.Correlation coefficients for maximum NDVI in a given year (year 0) with lagged 3-month running mean temperature (areaweighted and NPP-weighted).Positive correlations mean greater NDVI during (or following) warmer temperature.Filled circles indicate significance greater than the 95 % level.Shaded bars indicate the summer months (May-August).

Table 1 .
CO 2 observation stations included in each inversion.

Table 2 .
Trend and significance statistics for time series of interest, from 1986 through 2006."Trend" is the slope from linear least squares (LSQ) and Mann-Kendall (M-K) Sen slope; likewise, "Sig (p value)" is the p value from LSQ and M-K tests.Arctic zone and boreal zone are for BNA + BA; EU: Europe, NO: northern ocean.Italic values indicate 90 % significance level (p value 0.1 in at least one statistical test).

Table 3 .
Same as Table 2 but for the period from 1985 through 2012.