Reconciliation of top-down and bottom-up CO2 fluxes in Siberian larch forest

Carbon dioxide (CO2) fluxes by different methods vary largely at global, regional and local scales. The net CO2 fluxes by three bottom-up methods (tower observation (TWR), biogeochemical models (GTM), and a data-driven model (SVR)), and an ensemble of atmospheric inversions (top-down method, INV) are compared in Yakutsk, Siberia for 2004–2013. The region is characterized by highly homogeneous larch forest on a flat terrain. The ecosystem around Yakutsk shows a net sink of CO2 by all the methods (means during 2004–2007 were 10.9 g C m−2 month−1 by TWR, 4.28 g C m−2 month−1 by GTM, 5.62 g C m−2 month−1 and 0.863 g C m−2 month−1 by SVR at two different scales, and 4.89 g C m−2 month−1 by INV). Absorption in summer (June–August) was smaller by three bottom-up methods (ranged from 88.1 to 191.8 g C m−2 month−1) than the top-down method (223.6 g C m−2 month−1). Thus the peak-to-trough amplitude of the seasonal cycle is greater for the inverse models than bottom-up methods. The monthly-mean seasonal cycles agree among the four methods within the range of inter-model variations. The interannual variability estimated by an ensemble of inverse models and a site-scale data-driven model (the max-min range was 35.8 g C m−2 month−1and 34.2 g C m−2 month−1) is more similar to that of the tower observation (42.4 g C m−2 month−1) than those by the biogeochemical models and the large-scale data-driven model (9.5 g C m−2 month−1 and 1.45 g C m−2 month−1). The inverse models and tower observations captured a reduction in CO2 uptake after 2008 due to unusual waterlogging.


Introduction
Northern high-latitude ecosystems have experienced substantial perturbations due to the recent warming. Mean air temperatures increased by approximately 1.0 • C in the 50 • -60 • N region, which was twice the global mean increase of 0.5 • C during 1950-2015(Trenberth et al 2007, Hansen et al 2006. Climate models projected further increase in the Arctic temperatures up to 5 • C-7 • C by the end of this century (Collins et al 2013). The atmospheric CO 2 is increasing at a rate of about 2.0 ppm per year since year about 2001 (Patra et al 2005, Dlugokencky andTans 2017). Those changes will influence the carbon source and sink processes in the northern terrestrial ecosystems, however the trends and variabilities in net fluxes are highly uncertain at the moment due to the lack of observations. The high-latitude ecosystems store large amount of carbon as the soil organic matter in permafrost ( . It is also suggested that the enhanced CO 2 uptake in summer is potentially cancelled out by the release in autumn (Piao et al 2008).
Terrestrial biogeochemical models and atmospheric CO 2 inversions suggested no significant change in the annually-averaged carbon accumulation in the boreal forests in the past 20 years (Welp et al 2016). While the warming of the climate and rising atmospheric CO 2 concentration stimulates carbon uptake in summer, the soil carbon could be released in the form of CO 2 , leading to a positive climate-carbon feedback (Cox et al 2000, Koven et al 2011, Crowther et al 2016. For detecting the changes in net CO 2 balance in the boreal forests, measurements have been conducted for CO 2 fluxes at a larch forest at Yakutsk, Siberia using the eddy covariance method since early 2000s (Ohta et al 2014). Yakutsk is located in the center of east Siberia and the site is covered with the dominant vegetation type of deciduous needle-leaf (larch) forest that is one of the most vulnerable to climate change. This site is well suited for recording long-term changes in forest ecosystem uptake due to climate variations and change. Recognizing the importance of the region in global climate studies, measurements of atmospheric CO 2 concentrations in eastern Siberia are also ongoing using a network of tower sites (Japan-Russia stations; Sasakawa et al (2013)) and chartered aircrafts (Machida et al 2008).
Here we present an analysis of CO 2 fluxes from the bottom-up approaches and the top-down approach. The former approach includes the eddycovariance method (Ohta et al 2014), terrestrial ecosystem models from the GRENE-TEA (GRENE Terrestrial Ecosystem of the Arctic: the terrestrial research sub-project of the Arctic Climate Change Research Project of the GRENE (Green Network of Excellence) program by the Ministry of Education Culture, Sports, Science and Technology, Japan) model inter-comparison project (Miyazaki et al 2015), and a data-driven method using a machine learning algorithm (Ichii et al 2017). The top-down approach is the inverse modeling of the atmospheric CO 2 concentrations using atmospheric transport models. This study aims at direct comparison of top-down and bottom-up methods for clarifying the current status of CO 2 balance of larch forests in Siberia.

Data and methods
As a case study of the direct inter-comparison of the top-down and bottom-up methods, we selected an observation site at Yakutsk in Siberia, because of the availability of long-term CO 2 flux tower measurements and the highly homogeneous land surface condition, covered mostly by larch forests with almost flat topography (see figure S1 available at stacks.iop.org/ ERL/12/125012/mmedia, Friedl et al (2010)).
We used net CO 2 fluxes for the direct comparison, which denote net ecosystem production (NEP) for the bottom-up methods and net biome production (NBP, that is NEP plus disturbances) for the top-down method. The bottom-up methods are (1) the eddy covariance method using the tower observation (hereafter TWR), (2) the terrestrial biogeochemical models of GTMIP (hereafter GTM), and (3) the data-driven estimation using the support vector regression (SVR). The top-down method is the atmospheric inversion models (hereafter INV).
Areal means of 5 • × 5 • were used for SVR and INV for the comparison, because this spatial scale was the smallest scale in practical use of INV results (the scale issue is addressed later in detail). The monthly CO 2 fluxes were analyzed for seasonal and inter-annual variations.

Eddy covariance method
Tower observation of surface CO 2 flux was conducted at larch forest of Spasskaya Pad station (62.2550 • N, 129.2414 • E, 220 m amsl) near Yakutsk in eastern Siberia. CO 2 fluxes were measured at a 32 m tower site using the eddy covariance technique and archived at half-hourly time step through quality control and  That agrees with the general spatial scale of flux tower mesurements over forest, 100-1000 m depending on the atmospheric and surface conditions (e.g. Göckede et al 2004). Since the land cover around Yakutsk is fairly homogenous with an aerial coverage of about 500 km, in which 66% is deciduous needle-leaf forest (figure S1), we expect the tower flux measurements to represent a wider region.
The flux observation was limited to the warm season from May-September, and CO 2 fluxes during the other months were assumed only emission and estimated using temperature function (e.g. Lloyd and Taylor 1994). This study uses observed flux dataset for the period of 2004-2013.

Terrestrial biogeochemical models
Outputs from eight terrestrial biogeochemical models that participated in the GTMIP Stage 1 were used in this study (table 1; from Miyazaki et al (2015)). All the GTMIP models were driven by a common meteorological data. The model input data for the Yakutsk site were created from the nearest-point 6 hourly ERA Interim reanalysis data, which was then statistically fitted to the site using the 2005-2011 tower observation values (see table 2 in Sueyoshi et al (2016)). Models ran for the period from 1980-2013.

Empirical estimation: support vector regression
The empirical upscaling techniques provide another independent estimate of terrestrial CO 2 fluxes. We used the terrestrial net ecosystem exchange (NEE) empirically upscaled at 0.25 • spatial resolution and 8 day temporal resolution across Asia (Ichii et al 2017). The data used 54 eddy-covariance observation sites (278 site-years) and remote sensing products (e.g. land surface temperature, shortwave radiation, and water index) to estimate NEE based on SVR algorithm. 278 site-year data across Asia were used to train the model, and data from the Yakutsk Spasskaya Pad station were also included in the training process.
We used three different spatial coverage outputs of SVR-based CO 2 fluxes. One is based on 3 km × 3 km area centered on the eddy-covariance sites as site-level products using MODIS ascii subset as input (site-SVR). The other two estimates are based on spatial estimation of 0.25 • × 0.25 • spatial resolutions for continental distribution (continetal-SVR). We picked up 0.25 • × 0.25 • grid located at site, and 5 • × 5 • centered on the sites.

Atmospheric inversion estimate
We used CO 2 fluxes from four different inversion systems; global Eulerian-Lagrangian coupled atmospheric model (GELCA) based inversion system (GELCA-EOF; Zhuravlev et al (2013)), an atmospheric general circulation model-based chemistry transport model (ACTM) based inversion system (Saeki and Patra 2017), the nonhydorstatic icosahedral atmospheric model (NICAM)-based transport model (NICAM-TM) based inversion system (Niwa et al 2012), NIES's off-line atmospheric transport model (NIES-TM) based inversion system (Saeki et al 2013) (table 1). NIES-TM inversion system provides three different CO 2 fluxes using different observation network and inverse methodologies. Each inversion systems use different transport models, a priori surface flux dataset, atmospheric CO 2 observations (table S1), optimization method, and background errors, which causes variations in estimated a posteriori CO 2 fluxes. As mentioned above, obtained results from the inversions are monthly NBP, where the land biosphere fluxes include natural, land-use change, damages by insects and fire fluxes etc.

Mean annual balance and seasonal cycle
The mean annual, seasonal and monthly CO 2 fluxes during 2004-2007 were shown for each method in figure 1 when the results of all the models and observation were available. The annual fluxes were positive, i.e. carbon sink, for all methods (0.0109 for TWR, 0.00428 with the range from 0.000732 to 0.0119 for GTM, 0.00562 for SVR (3 km), 0.000863 for SVR, and 0.00489 with the range form −0.00275 to 0.00975 for INV; unit is kg C m −2 month −1 ), agreed among the methods within the ranges of GTM and INV models (shown by error bars in figure 1). The negative annual balance appeared in only one INV model, INV-4 (see figure S4(b) for detail). The seasonal fluxes also showed similar uptake in summer (June-August, JJA) and emissions in winter (December-February, DJF). As mentioned, the TWR measurements cover only for the period May-September. We estimated the cumulative CO 2 flux by TWR for these months alone, as 0.178 kg C m −2 (i.e. sink on land) while that for the extended period (October-April) was −0.0468 kg C m −2 (i.e. source to the atmosphere). Most of it (about 80%) was emitted in October and April when air temperatures were warmer than −10 • C. The fluxes in no-observation months accounted for about 70% of the annual cumulative flux (−0.0654 kg C m −2 during September-April), which was one third of the annual cumulative flux (+0.197 kg C m −2 during May-August). Our estimations of the fluxes in winter were similar to Wang et al (2011), in which the ranges of winter respiration in temperate and boreal ecosystems were shown to be about 40%-60% of the mean values. Since the contribution of winter respiration to the annual balance has been rather large and its uncertainty has been not negligible, further study on winter respiration is expected.
The mean absorption fluxes in June, July and August by TWR were 0.0853, 0.0740 and 0.0325 kg C m −2 month −1 . The ensemble-mean summer (June-August, hereafter JJA) absorption of INV (0.2236 kg C m −2 in JJA) was larger than those of GTM (0.0881 kg C m −2 in JJA), SVR (0.1019 kg C m −2 in JJA) and TWR (0.1918 kg C m −2 in JJA), indicating that the seasonal amplitudes by top-down models were larger than those of bottom-up models though its reason was remained open (see also figure S2). In contrast to the consistent CO 2 absorption in summer among all the models and observation, the variations of GTM in May and September-November are found to be both negative and positive. It suggests difficulty in representing a balance between photosynthesis and respiration in the seasons of transition between CO 2 uptake and release. Another possible reason is the simple representation of the GTM models considering either upper canopy or vegetation at the forest floor, only, though the structure of forest ecosystem is actually complex with various vegetation types.

Interannual variations in july
Since the observed CO 2 fluxes in July were most reliable (lowest gap in TWR data) during the absorption months, its interannual changes were examined. The CO 2 absorption flux was reduced after 2008 by TWR ( figure 2(a) The reduction of CO 2 fluxes after 2008 was not shown by the ensemble mean of GTM ( figure 2(a)). That is because the processes of vegetation damage caused by waterlogging have not been explicitly implemented. Nevertheless, two of the GTM models (GTM-4 and 5) have shown reductions after 2008 (table S2). Those changes are possibly induced by indirect effects of the unusual waterlogging, such as low soil temperature induced by excessive soil moisture. It should also be noted that GTM-1 and GTM-4 showed net release of CO 2 (i.e. negative values) in July fluxes for some years, which is inconsistent with observation showing net uptake during the summer months. In the case of GTM-1, the seasonal variation showed peak in CO 2 uptake in May-June but it decreased in July ( figure  S4(a)) because of high temperature causing increased respiration, resulting in low (and sometimes negative) NEP in July. Even with the net release of CO 2 in July, GTM-1 and GTM-4 produced net annual sink during 2004-2013. Another distinct point for GTM was the small interannual variation with its ensemble mean compared to TWR (figure 2(a)). It was clearly because the interannual changes of GTM models were out of phase each other. The divergence among the GTM models implies that important processes for interannual changes have been insufficiently considered.
On the other hand, the ensemble mean of INV showed similar interannual variation to that of TWR (figure 2(b), see table S3 for the averages and maximum-minimum ranges during 2004-2011). That is because the phases of interannual change have been roughly consistent among each INV model. The detailed reason for the similarities among the INV models have remained to be examined, but we confirm that some of the INV models employ similar inversion settings. The INV models used two sets of prior biosphere fluxes (CASA and VISIT), three different transport models (NIES-TM, NICAM-TM and ACTM) and three different inversion methods (GELCA-EOF, Kalman smoother, Matrix-based) (table S1). Moreover, INV-1 and INV-2 were looked to detect the reduction in CO 2 absorption after 2008 although the models did not use any observations from the Siberian region. The result suggested that the fluxes by INV around Yakutsk were constrained by the gas concentration data of remote observations, and might also that the reduction of CO 2 absorption after 2008 occurred on a greater regional scale.
The interannual variations simulated by SVR were smaller at any scale than TWR (figure 2(c)), particularly for continental-SVR (table S3). Besides, the reduction of CO 2 uptake after 2008 did not appear by any of SVR (figure 2(c)). Those disagreements were attributed, in part, to the tuning process of the SVR model objected to reproduce the 8 day satellite vegetation index. Namely, the important variables are chosen as explanatory variables to represent the 8 day vegetation index (e.g. land surface temperature, shortwave radiation, and water index), and the model structure (e.g. kernel parameters) are optimized to the 8 day vegetation index (Ichii et al 2017). As a result, the important factors or relations for interannual variation might not be included in the explanatory variables and in the model structure, hence leading to the weak sensitivity to interannual variations. Another possible cause is that the data period for tuning the SVR model is until 2008, leading to low representation of the reduction after 2008. The poor representation of the reduction in CO 2 absorption after 2008 by SVR apparently contradicts the implication from INV. However, if SVR model could present realistic variations and distributions after being tuned to interannual variations, much detailed discussion would be possible. It is thus needed to develop a tuning procedure for interannual variation in the SVR approach, e.g. by exploring explanatory variables important to interannual variations, and by examining to use anomaly from climatological seasonal cycle. It would also be effective to improve the regression procedure that is specific to the northern regions.

Difference in accordance with scale
In the comparison of top-down and bottom-up results, consideration for the difference in representative scale between the two should be of concern. As a starting point, differences of the mean seasonal cycle and the interannual variation at various scales were examined. The seasonal amplitude of site-SVR was closer to TWR than continental-SVR, naturally (figure 1). The summer absorption of site-SVR was smaller than TWR (figure S3), and that of continental-SVR is far smaller than TWR. That is probably because forest stands around the tower site are taller and their biomass is larger around the station than the surrounding forest (Dolman et al 2012). Another possible reason is a method of data screening, that the threshold to remove spike errors of observed CO 2 flux in preparing the reference data for machine learning of SVR is lower than that for eddy covariance estimation. Secondly, the interannual variation of the July CO 2 fluxes by site-SVR were closer to TWR in most of the years than those by continental-SVR (figure 2(c)). The average and maximum-minimum range during 2004-2011 by site-SVR agreed best with TWR among the SVR results (table S3).
To examine the correspondence to TWR in accordance with scale, the monthly fluxes by SVR and one of the INV models (GELCA-EOF) averaged at multiple scales were compared with TWR. The scatter plots and regression lines for SVR ( figure 3(b)) showed that the monthly fluxes became larger and closer to TWR as the representing scale became small (from 10 • × 10 • to 3 km × 3 km) in the positive quadrant (i.e. sink). That is natural because the SVR model was established at 3 km scale using the observed data at the site and applied to a continental scale using the largescale datasets, and because the forest biomass and hence the CO 2 absorption around the tower station was rather large as mentioned above. Similarly to SVR, the monthly fluxes by GELCA-EOF ( figure 3(a)) averaged in 5 • × 5 • became larger than those averaged in 10 • ×10 • in the positive quadrant. However, the dif- ferences from TWR were larger for those averaged in 5 • ×5 • than in 10 • × 10 • . That could be attributed to the overestimations of CO 2 sink by GELCA_EOF in summer. There are a few possible reasons for the overestimation, e.g. (1) results of an inversion model generally include numerical fluctuation at a small scale resulting in large values than actual at a small scale, (2) the averages for 10 • × 10 • could pick up the grids with less variation in a priori field than those for 5 • × 5 • ( figure S4(d)).
As shown from the comparison above, the estimations of CO 2 fluxes by the top-down method (one of INV models) and a bottom-up method (SVR) presented the consistent relation to tower observation in accordance with scale. Such comparison becomes possible because CO 2 fluxes are estimated under a common setup of the SVR model both at a site scale and at a large scale. The comparison was limited to a simple one in this study, since the object of the SVR model (to reproduce 8 day satellite vegetation index) was different from the target phenomena (seasonal cycle and interannual variations). Actually, the estimates by the SVR model and by the inversion model agree well at sub-continental scale in mid and high latitudinal regions (Kondo et al 2015). It hence implies that SVR has a potential to provide a comprehensive understanding to the scale issue. Detailed interpretations would be realized if the SVR model is optimized to a target phenomenon and analysis is made at a suitable scale for the counter-method of the comparison, e.g. comparisons between TWR (or GTM) and SVR at a site scale, and between INV and SVR at a continental scale. Furthermore, examinations including spatial distributions and temporal changes would be useful to understand the different representations between top-down and bottom-up methods.

Conclusions
This study compared estimates of CO 2 fluxes by tower observation, top-down inverse models and two bottom-up methods (biogeochemical models, and a data-driven model) in Yakutsk, Siberia for 2004-2013. The region is characterized by highly homogeneous larch forest on a flat terrain. The ecosystem around Yakutsk shows a net sink of CO 2 by all the methods (0.0109 by tower observation, 0.000863-0.00562 by the bottom-up models, and 0.00489 by the top-down method; units are kg C m −2 month −1 ). The mean seasonal cycles during 2004-2007 showed good agreement among the methods, within the range of variation of each model of inversion and biogeochemical models. Mean absorption in summer (June-August) ranged from 0.0881-0.1918 kg C m −2 month −1 by three bottom-up methods and was 0.2236 kg C m −2 month −1 by the topdown method, showing that the seasonal amplitude by top-down models was larger than that of bottom-up models.
The interannual variation in July by TWR presented that the CO 2 flux was reduced after 2008 because of damage of forest stands by the unusual waterlogging. However, similar changes did not appear by most of the GTM since such processes have not explicitly been considered. Besides, the interannual variance of the ensemble mean of GTM was much smaller than TWR. That was because the variation of each model was out of phase and canceled out each other. It implies the needs to examine important factors and processes for interannual variations by examining GPP and respiration, which will lead to improve the representations of carbon flows and allocations in GTM models. In contrast, the interannual variance of the ensemble mean of INV was similar to TWR, because the phase of variation by each INV model rather agreed with each other. The reason for the similarities among the INV models have remained to be examined. The reduced uptake after 2008 was presented by INV though the results of only two models were available. It suggested that fluxes by INV around Yakutsk are constrained by the gas concentrations from the remote observations. The interannual variations simulated by SVR were smaller than that observed (TWR) and the observed reduction after 2008 was not presented by SVR, since the SVR model has not been established for interannual variations. It will be needed, for improving the interannual representation, to include important variables for interannual variations as the explanatory variables and to optimize the model structure (e.g. kernel parameters) at the target time scales.
The comparison of TWR with site-based and largescale SVR revealed that the site-SVR corresponds better with TWR, showing the local influence of the exceedingly large forest around the tower site. The comparison of TWR with one of the INV models at two scales showed similar correspondence. If the SVR model is established for a targeted issue (e.g. interannual variation), its result could be examined in extended ways, e.g. to make a quantitative comparison, to discuss horizontal distribution, etc Furthermore, comparisons for other variables such as GPP and respiration and for different regions are desired to understand reasons for the scale gap and hence to reduce the uncertainty in terrestrial carbon balance.