Soil CO2 Emission Largely Dominates the Total Ecosystem CO2 Emission at Canadian Boreal Forest

The natural carbon dioxide (CO2) emission from the ecosystem, also termed as the ecosystem respiration (Reco), is the primary natural source of atmospheric CO2. The contemporary models rely on empirical functions to represent decomposition of litter with multiple soil carbon pools decaying at different rates in estimating Reco variations and its partitioning into autotrophic (Ra) (originating from plants) and heterotrophic (originating mostly from microorganisms) respiration (Rh) in relation to variation in temperature and soil water content. Microbially-mediated litter decomposition scheme representation are not very popular yet. However, microbial enzymatic processes play integral role in litter as well as soil organic matter (SOM) decomposition. Here we developed a mechanistic model comprising of multiple hydro-biogeochemical modules in the soil and water assessment tool (SWAT) code to explicitly incorporate microbial-enzymatic litter decomposition and decomposition of SOM for separately estimating regional-scale Ra, Rh and Reco. Modeled annual mean Reco values are found varying from 1,600 to 8,200 kg C ha−1 yr−1 in 2000–2013 within the boreal forest covered sub-basins of the Athabasca River Basin (ARB), Canada. While, for the 2000–2013 period, the annual mean Ra, Rh and soil CO2 emission (Rs) are varying within 800–6,000 kg C ha−1 yr−1, 700–4,200 kg C ha−1 yr−1 and 1,200–5,000 kg C ha−1 yr−1, respectively. Rs generally dominates Reco with nearly 60–90% contribution in most of the sub-basins in ARB. The model estimates corroborate well with the site-scale and satellite-based estimates reported at similar land use and climatic regions. Mechanistic modeling of Reco and its components are critical to understanding future climate change feedbacks and to help reduce uncertainties particularly in the boreal and subarctic regions that has huge soil carbon store.


INTRODUCTION
The natural carbon dioxide (CO 2 ) emission from the ecosystem, also termed as the ecosystem respiration (R eco ) is the primary natural source of CO 2 throughout the globe (Ciais et al., 2014). There is an urgent need to completely understand the soil carbon conversion processes and the associated land-surface gas exchange in estimating climate change feedback, which ultimately leads to the release of soil carbon as CO 2 to the atmosphere (Mitchard, 2018). In natural ecosystems, R eco can be partitioned into the soil heterotrophic respiration (R h ) that mostly resulting from soil microorganisms and the autotrophic respiration (R a ) resulting from plant species (both from aboveground components and roots) (Hicks Pries et al., 2013;Hicks Pries et al., 2016). Globally, R a is the predominant component of the terrestrial carbon budget with photosynthetic carbon consumption rate of 54-71% (Ryan et al., 1997). In permafrost regions, R a accounts for 40-70% of the total ecosystem respiration (Hicks Pries et al., 2013).
The R a and R h have different adaptation mechanisms in response to climate change. For example, rising temperature can accelerate microbial enzymatic activities and thus enhance R h and facilitate nutrient dynamics based on substrate availability (Davidson and Janssens, 2006;Manzoni et al., 2008). On the other hand, rising temperature has different influence in plant growth. The aboveground component on R a shows positive response while root shows near neutral response (Hick Pries et al., 2013). For example, 2.3°C warming of soil leads to 20% increase in aboveground productivity at an Alaskan site (Natali et al., 2012). With an increase of 2°C, root respiration didn't show any significant changes, while heterotrophic respiration did 21% increase linked with increased microbial activities (Wang et al., 2014). Therefore, quantifying the distinct contributions of R eco , R a , and R h is necessary to estimate climate feedbacks and sensitivity in the rapidly changing subarctic regions but state of the art monitoring techniques, such as eddy covariance and remote sensing techniques, are still unable to directly partition R a and R h from R eco (Davidson and Janssens, 2006;Konings et al., 2019). In addition to the distinct characteristics and feedbacks to climate change processes by R a and R h , respiration partitioning is particularly important as it has been reported that global land carbon sink has been increasing in recent years (Ciais et al., 2019).
Climate change-linked future projections of ecosystem respiration are challenging because estimating litter and soil organic matter (SOM) decomposition rates and anticipated soil derived CO 2 feedbacks resulting from anthropogenic warming are both seen as being highly uncertain (Collins et al., 2013;Crowther et al., 2016). Many experimental and modelling studies have been performed to determine the R eco and its partitioning and gross primary production (GPP) at different land use types (Hardie et al., 2009;Hicks Pries et al., 2013;Senapati et al., 2018). On the other hand, predicting the net ecosystem exchange (NEE), and R eco require the development of sophisticated land-surface models. These modeling approaches can be categorized into: 1) agroecosystem models (models used to simulate agricultural system functioning), and 2) the Earth System Models (ESMs -that are used to simulate Earth system processes). Both modelling approaches commonly use empirical formula of multiple soil C pools decaying at different rates to calculate R eco (Del Grosso et al., 2005;Davison and Janssens, 2006;Oleson et al., 2010;Clark et al., 2011). The CENTURY model used simplified functions of soil temperature and moisture for estimating soil respiration and its partition to autotrophic and heterotrophic components (Del Grosso et al., 2005). The decaying module of multiple soil C pools in Daily CENTURY model (DayCent) has also been incorporated into Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) (Qiu et al., 2018) and Community Land Model (CLMs) (Lawrence et al., 2019). The microbial processes in ESMs are simplified into linear, empirical equations (Crowther et al., 2014;Crowther et al., 2019). It is found that contemporary ESMs cannot reproduce grid-scale variation in soil C due to missing key processes and the predicted global carbon stocks in the fifth Coupled Model Intercomparison Project (CMIP), leading to 6fold difference in predicted data (Todd-Brown et al., 2013). Microbial activity responses (as reflected in soil heterotrophic respiration) are expected to increase with warming (Karhu et al., 2014;Walker et al., 2018). However, the relationship between warming and soil carbon loss, as well as overall ecosystem respiration, is not straightforward as the microbially mediated litter and SOM decomposition processes are not linearly correlated with temperature (Melillo et al., 2017). Additionally, microbially-mediated decomposition of SOM is not only an important biological-driven process for carbon conversion but also play a major role in overall global nutrient dynamics (Manzoni et al., 2008).
Apart from the environmental variables, the turnover of organic material is directly controlled by soil microbes (Crowther et al., 2019). As a result, there is a need to improve microbial processes for global carbon modeling estimates (Wieder et al., 2013;Crowther et al., 2019). Microbial abundances are relatively larger in arctic and subarctic regions (Serna-Chavez et al., 2013;Xu et al., 2013). They are also responsible for biogeochemical cycling of nutrients (Crowther et al., 2019). However, different types of microbes require favorable soil redox conditions for their growth (DeAngelis et al., 2010). Thus, soil redox condition is an essential measure for the dynamics of nutrients as well as soil organic matter, which microbes used as a substrate (Bhanja et al., 2019a;Bhanja et al., 2019b). It has been found that nutrient availability directly controls soil organic matter stock and terrestrial carbon sink (Wieder et al., 2015). Alteration of hydrological processes has profound impact in soil carbon mineralization (Anthony et al., 2018). Therefore, integration of hydrological processes along with redox/biogeochemical and microbial processes would definitely improve model estimates (Bhanja et al., 2019a;Bhanja et al., 2019b). Global-scale respiration studies continue to be sought for further improving the modeling processes and soil respiration estimates (Todd-Brown et al., 2013;Luo et al., 2016;Wang et al., 2020;Wang et al., 2021). A decomposition feedback to warming requires accounting explicitly for not only temperature but also nutrient availability and microbial activities (Davidson et al., 2012).
Substantial amount of soil organic carbon can be released due to enhanced microbial activities associated with climate warming in arctic, subarctic regions (Schuur et al., 2015). The arctic and subarctic regions are the most sensitive regions to the combined climate change impacts due to declining permafrost, glacial retreat and change in freeze-thaw cycles on its ecosystems (Bates et al., 2008). An abrupt thawing in lakes due to global warming can lead to faster mobilization of deeper stored carbon (Anthony et al., 2018). Therefore, in these regions, modeling the carbon mobilization processes and their feedbacks are much sought after (Schuur et al., 2009). The present study describes the development of a mechanistic model to separately simulate regional-scale, autotrophic, heterotrophic and total ecosystem respirations. This being achieved through an explicitly integrated inclusion of microbially-mediated decomposition of litter as well as SOM, redox processes and hydrological processes. The model performance was assessed at the boreal forest covered Athabasca River Basin (ARB), Canada (Fig. S1) using both sitescale, as well as remote-sensing data.

Model Description
The Soil and Water Assessment Tool (SWAT) is widely used for simulating regional-scale hydrology due to its comparatively lower computational requirements (Arnold et al., 1998). Integrated hydrobiogeochemical modeling was performed at a daily time step. Our microbial kinetics and thermodynamics (MKT) model is integrated within SWAT framework, known as SWAT-MKT (Bhanja et al., 2019a;Bhanja et al., 2019b;Bhanja and Wang, 2020). While the SWAT calculates hydrological processes, MKT evaluates microbially-mediated dynamics of biogeochemical processes in soils. At daily time step, two models exchange data and update their variables. In this process, MKT rely on daily soil moisture and temperature simulation from SWAT and independently simulate its biogeochemical processes. The oxygen diffusion and related oxidizing reactions are also modelled. Oxygen diffusion processes are described in supplementary information.
In order to simulate the processes of multiple chemical reactions, microbial kinetics and thermodynamics are simultaneously considered in SWAT-MKT model. The energy balance and sequence of the chemical reaction to be simulated is modeled through thermodynamics approach following the Nernst equation. Considering a chemical reaction: Where, A, Aare the normal and reduced form of the chemical component modeled, respectively. D and D + are the chemical component that donate electron in the reaction and its oxidized form, respectively. x, y, p and q, are the stoichiometric coefficients of reactants and products.
The redox potential (E h ) dynamics was modeled following the Nernst's equation (Stumm and Morgan, 1996).
Where, E h and E 0 are the redox potential (mV) and the standard electrode potential (mV), respectively. R is the universal gas constant; T is the ambient temperature in Kelvin; n is the number of electrons transferred in the chemical reactions; F and Q are the Faraday's constant and the reaction quotient, respectively. '[ ]' represents concentration in molality and 'γ' is the activity coefficients in molality. The chemical reactions in soil-water zone occur at a sub-daily time step catalyzed by various microbes. Otherwise, if the reactions would have happen in thermodynamic control only, the rates would be much slower e.g. in microbe limiting conditions, in groundwater systems. The microbes catalyzed the reactions after taking energy from the system. We have considered this non equilibrium condition by after computing and removing the microbial energy requirement from the system and updated the net redox potential at a daily time step (for details, please check Bhanja et al., 2019a).
The mass balance of different chemical species is achieved through modeling microbial kinetics for estimating the reaction rates of individual reactions at a daily time step. Dual Michaelis-Menten kinetics is followed here, the reaction rate can be represented as: .
Where, rate of the reaction is represented by R in molal/day; specific microbial activity as Q max in mol/mol of biomass per day; B is the microbial biomass in mole biomass/l of water; [D] is the concentration of the electron donor species in Eq (1) in molal; [A] is the concentration of the electron acceptor species in molal, and K D and K A are the half saturation constants for the electron donor and acceptor species, respectively in mole/l. K D and K A are otherwise termed as the Michaelis-Menten constants.
Major soil oxidation-reduction reactions considered in this approach were given in Supplementary Table S1 and their reaction quotient values are provided in Supplementary Table  S2. The reactions are shown in Supplementary Figure S3 and further details are provided in Bhanja et al. (2019a); Bhanja et al. (2019b) and Bhanja and Wang (2020).
New carbon cycle capabilities are incorporated into SWAT to simulate ecosystem respiration components from litter decomposition, root respiration, above ground respiration and respiration component from the dissolved organic carbon transformation by enzymatic processes (Supplementary Figures S2, 3). The entire chemical processes considered in the new version of SWAT-MKT are shown in Supplementary Figure S3.
Various spatial and meteorological datasets are used to builtup the SWAT model for the ARB. Shuttle Radar Topography Mission (SRTM)-based Digital Elevation Model (DEM) data, at a resolution of 90 m × 90 m, are obtained from the Consultative Group on International Agricultural Research (CGIAR) (Jarvis et al., 2008). Land-use data are obtained at 1 km × 1 km, from the International Geosphere-Biosphere Programme, Data and Information Systems (IGBP-DIS) initiative archived at USGS (Loveland et al., 2000). The soil map is obtained from the Agriculture and Agri-Food Canada at 1:1 million resolution (SLC, 2010). Elevation ranges from 207 to 3,669 m in the ARB. We have defined 320 soil classes and 11 land-use classes. Watershed delineation through SWAT with a 200 km 2 threshold gives rise to 131 sub-basins within the ARB. Hydrologic Response Units (HRU) are defined using four slope classes: 5, 10, 15 and 20%, and with 10, 5 and 10% thresholds for land-use, soil and slope, respectively. A total of 1,370 HRUs are obtained through this process in this the region. Meteorological data are used at a daily scale for precipitation, minimum and maximum temperature, were obtained from the GoC (2016) database. Solar radiation, relative humidity and wind speed data are obtained for 230 stations from CFSR (CFSR, 2016). Details on these data processing and model built-up can be found in Shrestha et al. (2017).

Soil Mineralization and Litter Decomposition
The equations should be inserted in editable format from the equation editor.
Soil mineralization module is modeled considering two carbon pools: active and passive. The active pool represents fraction of active litter including microbial biomass (Fujita et al., 2014). We used microbial enzymatic litter transformation approach that is an advancement of the CENTURY model's simple, first order kinetics-based litter decomposition approach (Parton et al., 1987;Parton et al., 1994;Fujita et al., 2014): Where, litter decomposition rate originally adopted in CENTURY model: R d i,C (gC kg −1 soil d −1 ). i represents the substrate type i.e. active and passive. C i (gC kg −1 soil) is the carbon content within active or passive substrate.
Litter decomposition can also be modeled through microbial enzymatic approach following one-substrate Michaelis-Menten kinetics (Fujita et al., 2014). Thus, the new decomposition rate becomes: Where, Km i is the half-saturation constant or Michaelis-Menten constant. k i,M is the decomposition coefficient of C i and can be estimated separately for the active (AC) and passive (PA) substrates as: Where, k AC,C and k PA,C are decomposition coefficients used in CENTURY model for the active and passive substrates, respectively. Km AC is the half-saturation constant for active substrate and is approximated as 0.3 g C kg −1 soil (Allison et al., 2010). Km PA is the half-saturation constant for the passive substrate and its value is taken as 600 g C kg −1 soil (Allison et al., 2010). C b (g C kg −1 soil) represents microbial biomass and its value is approximated as the median microbial biomass (0.87 g C kg −1 soil) from a global-scale study of Cleveland and Liptzin (2007). C T represents total carbon stock of soil and its value is approximated as the global total soil carbon (46 g C kg −1 soil; Cleveland and Liptzin, 2007). If considering microbial biomass being an active component of the litter, microbial biomass decomposition can directly be assumed to be proportional to the litter decomposition. The new litter decomposition rate (R d i,MM ) becomes (Fujita et al., 2014): Actual soil respiration rates (R LD ) from litter transformation is estimated as follows (Fujita et al., 2014).
Where, e i,m represents the growth efficiency of microbes during assimilation of either active or passive substrates and its value is estimated as 0.45 (Fujita et al., 2014). I m,c is an decomposition inhibition factor (its value varies from 0 for full to one for no inhibition) and at present its value taken as one also resembles the CENTURY model parameterization (Fujita et al., 2014). O m,c is the overflow of carbon due to limiting nitrogen concentration and its value is taken as 0 without proper data to represent the processes.

Root Respiraton
Root respiration (R r ) is an essential component of soil respiration. However, SWAT does not have the ability to simulate root respiration. To simulate root respiration, we have developed a new sub-module within SWAT following Li et al. (1994): Where, CO 2 produced by roots for nitrogen uptake: R n (13.8 mg C meq −1 N; Veen, 1981;Li et al., 1994). Nitrogen uptake rates of plant is represented as U n (kg N ha −1 d −1 ). CO 2 produced by roots due to their growth: R rg (19.19 mg C g −1 dry matter; Veen, 1981;Li et al., 1994). Root biomass growth at a day: BG r (g dry matter ha −1 ). CO 2 produced as a function of root maintenance: R rm (0.288 mg C g −1 dry matter d −1 ; Veen, 1981;Li et al., 1994). B lr is the living root biomass (g dry matter ha −1 ).

Satellite-Based Estimates of Autotrophic Respiration
We used gross primary production (GPP) and net primary productivity (NPP) data from the observation of the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors (Zhao et al., 2005;Zhao et al., 2006;Zhao and Running, 2010). Annual mean MOD17 products are used at a spatial resolution of 30 arcsec. MOD17 is the first satellite derived continuous data Frontiers in Environmental Science | www.frontiersin.org June 2022 | Volume 10 | Article 898199 product for vegetation productivity at the global-scale . The algorithm includes several satellite derived parameters such as the land cover, fractional photosynthetically active radiation and leaf area index along with the meteorological variables . NCEP/ DOE reanalysis II outputs are used for meteorological parameters (Zhao and Running, 2010). Detailed descriptions of the MOD17 products can be found in Zhao et al. (2005); Zhao et al. (2006). Satellite-based R a is estimated by subtracting NPP from the GPP data (Bond-Lamberty et al., 2018).

Assumptions and Limitations
In order to compute the total soil respiration, we only used the heterotrophic and autotrophic components. Here, we have not considered the geological CO 2 emission (Andrews and Schlesinger, 2001). Ecosystem respiration also includes animal respiration. However, due to the cold climatic conditions and presence of very low number of animals at the ARB (Weber et al., 2015), the respiration from animals are not considered at present in our approach. Other assumptions and limitations associated with the basic version of the model can be obtained from Bhanja et al. (2019a); Bhanja et al. (2019b) and Bhanja and Wang (2020).
In natural conditions, all of the soil microbes do not produce enzymes or produce at a slower rate to take part in the decomposition activities. These types of microbes restrict/slow down the decomposition process (Kaiser et al., 2015).

Heterotrophic, Autotrophic and Total Ecosystem Respiration
Annual mean ecosystem respiration (R eco ) did show spatial variability, however, most of the predicted values for the 2000-2013 period are in the range of 1,600-8,200 kg C ha −1 yr −1 (Figure 1). The R eco estimates were compared with the available site-scale measurements from the Canadian boreal forest locations data retrieved from FLUXNET 2015 (Pastorello et al., 2017). Most of the site-scale R eco values vary within 4,000-10,000 kg C ha −1 yr −1 , but with values on three sites exceeding 11,000 kg C ha −1 yr −1 (Supplementary Table S3). The modeled R eco estimates were in lower range compared to the other boreal forest R eco observations shown in Supplementary Table S3. These reported sites are however located either at more southern or the same latitude as our study area. Therefore, overall the prevailing climatic conditions at the ARB are more nudging towards a arctic-like climate, with more limited respiration rates. In general, the respiration values were lower during winter months, both climatic factors and the presence of deciduous trees do account for these lower rates (Cumming, 2001). The mean heterotrophic respiration (R h ) values did mostly vary from 700 to 4,200 kg C ha −1 yr −1 in 2000-2013 at the ARB (Figure 1). The data showed strong  Table S4) retrieved from the SRDB archives (Bond-Lamberty and Thomson, 2010). The Quebec sites are located far southern areas with little warmer climate and thus exhibiting higher R h . Most of these site-scale values varied between 1700 and 5,900 kg C ha −1 yr −1 , but on a few occasions (Quebec sites) with values of >10,000 kg C ha −1 yr −1 were recorded (Supplementary Table S4). Mean root respiration (R r ) values (300-2,800 kg C ha −1 yr −1 ) were however found to be lower than the mean R h values in our study at ARB (Figure 1). In general, R h is found to be higher than R r in global boreal sites (Bond-Lamberty and Thomson, 2010) and at an Alaskan site (Hicks Pries et al., 2013). R h is not only contributing higher toward R s , the contribution rate has been  increased from 54 to 63% over the years during 1990-2014 at a global study (Bond-Lamberty et al., 2018). Annual mean autotrophic respiration from above-ground vegetation components (R abv ) did vary from 200 to 2,900 kg C ha −1 yr −1 . In general, the combination of R r and R abv that is the mean autotrophic respiration (R a ) is estimated to be higher than the mean R h (Figure 1). Annual mean R a varies within 800-6,000 kg C ha −1 yr −1 (Figure 1). Most of the ARB is covered by forests (Bhanja et al., 2018), which may account for the higher magnitude of autotrophic respiration compared to its heterotrophic counterpart. Contribution of above-ground canopy respiration to R eco is significant at boreal forests (Ryan et al., 1997) and at arctic climates (Hick Pries et al., 2013). Thus R abv along with the R r , makes R a dominant over R h in forest. The R a dominancy was also observed in savanna and grasslands (Ma et al., 2007), and in peatland ecosystems (Hardie et al., 2009). The autotrophic and heterotrophic respiration estimates are well in line with the data reported in Goulden et al. (2011) from Canadian boreal forests with values reported for R a was estimated at 2000-4,500 kg C ha −1 yr −1 and R h at~2000 kg C ha −1 yr −1 . R a reported by Ryan et al. (1997) at eight Canadian boreal forest sites with values between 3,120 and 6,110 kg C ha −1 yr −1 are also matching our modelled estimates. Bond-Lamberty et al. (2010) reported R h within 200-6,000 kg C ha −1 yr −1 at boreal locations and 100-900 kg C ha −1 yr −1 at arctic locations across the globe. In general, mean soil respiration (R s ) also show spatial patterns with values from 1,200 to 5,000 kg C ha −1 yr −1 at the ARB (Figure 1).
The lower values of R s were mainly found at the Southern ARB subbasins dominated by mountains. The R s values do also well comparable with the site-scale estimates from many of the other Canadian boreal sites except those at Quebec (located in southern latitude with warmer climate) and some Saskatchewan sites (Supplementary Table S4).

Relationships Between the Respiration Components
Relationship between R h and R s shows near equal contribution of soil autotrophic and heterotrophic respiration to total soil respiration in most of the sub-basins (Figure 2). Root contribution to total soil respiration (RC) values did varies from 0.1 to 0.6 (occasionally 0.7). This generally matches the field-scale estimates, which further increases our confidence in the modelled estimates (Supplementary Table S4). The ratio of R a to R eco did show varying contribution of R a to R eco from 30 to 80% in most of the sub-basins ( Figure 2). The estimates are well within the ranges reported in previous studies (40-80% in Nowinski et al., 2010;40-70% in Hicks Pries et al., 2013) at similar eco-climatic regions. R s contributes to nearly 60-90% of R eco in most parts of the study area ( Figure 2). Contribution of R r to R eco lies within 10-45%. Results are consistent with the observation of Hicks Pries et al. (2013) at arctic climate (15-35% contribution reported). ARB is mostly covered by forest (Supplementary Figure S4) and occurrence of comparatively lower annual mean soil temperature (<2°C, Supplementary Figure S5) are the two main reasons for the dominance of autotrophic respiration toward the total ecosystem respiration (Ryan et al., 1997;Hicks Pries et al., 2013;Crowther et al., 2016). The relative proportions of R h to R eco were also found to be consistent with the values reported in previous studies (Hardie et al., 2009;Schuur et al., 2009;Hicks Pries et al., 2013).
The respiration partitioning and their ratio show some interesting facts. Our work shows that the relationships between R s and R eco and R a and R eco follow linear relationship (r 2 > 0.74, p < 0.001, Figures 3B, C). Relationship between R h and R s follow linear pattern (r 2 = 0.61, p < 0.001, Figure 3A). Similar relationships are also observed by Bond-Lamberty et al. (2004) in different locations across the globe.
Global anthropogenic CO 2 emissions (combination of fossil fuel and land use change) has been increased from 4.5 to 11 Gt C yr −1 from 1960-1969-2018(Friedlingstein et al., 2019. Global terrestrial ecosystem carbon sink has increased from 1.3 to 3.2 Gt C yr −1 from 1960-1969 to 2009-2018 and subsequently slowing down the atmospheric CO 2 concentration increase (Friedlingstein et al., 2019). Terrestrial ecosystems are acting as a carbon sink for approximately 29% annual anthropogenic CO 2 emissions during the last decade (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and the magnitude is higher than the ocean sink rates (~23%) (Friedlingstein et al., 2019). Although the direct link is unclear, it has been reported that atmospheric CO 2 concentration is also sensitive to terrestrial water storage change at global-scale with declining values associated with rapid increase of CO 2 concentration (Humphrey et al., 2018).

Modelled Estimates Comparison With Remote Sensing Data
Satellite-based estimates of R a ( Figure 4A) show similar spatial patterns on comparing with the modeled R a (Figure 1). Satellitebased R a varies from 800 kg C ha −1 yr −1 to as high as 3,943 kg C ha −1 yr −1 across the sub-basins of ARB. In general, modeled R a is aligned with the satellite-based R a in the mid-R a region (1,500-4,000 kg C ha −1 yr −1 ; Figure 4B). The dissimilarity, in lower and higher R a ranges are result of various well known issues with the satellite-based approach. Several studies have reported erroneous satellite-based NPP estimates (~15% less estimates comparing the observations). These were found to be associated with the interference from the autotrophic respiration estimation from the neighboring areas (Ito, 2011). This make the satellite-based R a value smaller than its real value. The GPP and NPP database developed using different available meteorological datasets are also showing an overestimation of the indices when comparing with the GPP (~30% higher GPP reported using NCEP data) and NPP (15-20% higher NPP using NCEP data) developed using observed meteorological data . The NCEP/DOE reanalysis II meteorological data were used to develop the global-scale GPP and NPP products (Zhao and Running, 2010)-this can also be a further reason for the overestimation. Turner et al. (2006) reported overestimation of MODIS NPP and GPP products at regions with comparatively lower productivity e.g. Boreal forest regions.

CONCLUSION
Although the terrestrial ecosystem respiration is one of main components of climate change feedbacks, the sign and the magnitude of this feedback is highly uncertain in future. Our approach using a new integrated hydro-biogeochemical model show reasonable estimates of regional-scale R eco and its subsequent partitioning into R a and R h on comparing with the available data at the boreal forest covered Athabasca river basin, Canada. Annual mean R eco ranging between 1,600 and 8,200 kg C ha −1 yr −1 in 2000-2013. The R s was dominated and contributed 60-90% toward R eco . The model estimates are in line with the site-scale measurements reported at similar land use and climatic regions. Satellite-based estimates of R a also show similar patterns as of the modeled estimates.
Reasonable performance of our model for simulating regionalscale R eco and its components show that the bottom-up approach performs good for estimating the respiration components. Anthropogenic CO 2 emission is a magnitude lower than the global ecosystem respiration. The accurate quantification of respiration components warrants better understanding of their feedbacks associated with the anthropogenic warming. Separate estimation of the respiration components is getting more attention here due to their differential responses to warming. It would be interesting to integrate the microbial enzymatic kinetics-based approach with widely used ESMs as a core module of soil greenhouse gas emissions.

AUTHOR CONTRIBUTIONS
SB acquired the data, developed the model code and completed the analyses with inputs from JW and SB wrote the manuscript with inputs from JW and RB.

ACKNOWLEDGMENTS
This is a short text to acknowledge the contributions of specific colleagues, institutions, or agencies that aided the efforts of the authors. SB. acknowledges support from the Indian Institute of Science in the form of C. V. Raman Postdoctoral Fellowship for partly carrying out the study. JW. would like to thank the Alberta Economic Development and Trade for the Campus Alberta Innovates Program Research Chair (No. RCP-12-001-BCAIP). The authors acknowledge FLUXNET and SRDB network for making their data available to public. We appreciate the MODIS land team for the MOD17 GPP and NPP data.