Influence of environmental factors on autotrophic, soil and ecosystem respirations in Canadian boreal forest

Ecosystem respiration (R eco ) and its components, the autotrophic respiration (R a ) and soil respiration (R s ) are the essential indicators of the global carbon cycle. They are represented as functions of either temperature or soil moisture, or a combination of both in the widely-used Earth System Models (ESMs). Thus, it is difficult to evaluate the influence of other environmental factors (such as, precipitation, soil temperature, dissolved oxygen level and oxidation reduction potential (ORP)) on R a , R s and R eco . Here we introduced microbially mediated, detailed carbon cycle processes within our mechanistic model to address this issue. Dominance analysis using a multivariate approach was performed to find out the influence of individual environmental factors on R a , R s and R eco in the cold climate regions of Athabasca River Basin (ARB), Canada. Contribution of the 6 predictor vari- ables, including air temperature, precipitation, soil temperature, water-filled pore space (WFPS) used as a proxy of soil moisture, dissolved oxygen level, and ORP, on R a , R s and R eco were estimated based on the R 2 values originated from multiple regression analyses. Our results showed that the prevailing temperature (both air and soil) and dissolved oxygen levels are the major influencing factors on R a , R s and R eco . WFPS is found to be the least influential factor on respiration estimation. Output of this study can be used to consider the crucial roles of environmental drivers in R a , R s and R eco estimation in the development of future ESMs.


Introduction
Ecosystem respiration (R eco ) is primarily responsible for soil carbon loss and natural carbon dioxide (CO 2 ) emission (Ciais et al., 2014). Autotrophic respiration (R a ) from the vegetations and heterotrophic respiration (R h ) from soil microorganisms are majorly constituting R eco in natural ecosystems (Hicks Pries et al., 2013;Hicks Pries et al., 2015). Site-scale estimation of R eco (Chambers et al., 2004;Knohl and Buchmann, 2005;Hardie et al., 2009;Nowinski et al., 2010;Hicks Pries et al., 2013), R a (Chambers et al., 2004;Bond-Lamberty et al., 2004;Bond-Lamberty and Thomson, 2010b), R h (Chambers et al., 2004;Bond-Lamberty et al., 2004;Bond-Lamberty and Thomson, 2010b;Bond-Lamberty et al., 2018) are well studied across the globe. Relationship between R eco and soil temperature and moisture are reported in earlier studies (Knohl and Buchmann, 2005;Misson et al., 2007). R s variations with temperature and soil moisture Thomson, 2010a, 2010b;Hursh et al., 2017) were observed in different parts of the globe. However, relationships between R eco , R s and R a with the environmental (temperature, precipitation, soil temperature, WFPS) and chemical environmental drivers (dissolved oxygen level, ORP) are not readily available.
Estimating the influence of environmental factors on R eco is a necessary step to understand and manage future CO 2 emissions particularly at the cold regions, where soil freeze and thaw cycle change associated with global warming could cause a shift in carbon mobilization mechanism in the near future (Bhanja and Wang, in-Press). The ongoing carbon cycle processes are changing in cold regions due to climate change linked decline of permafrost, glacial thinning and the changing pattern of freeze thaw cycle (Bates et al., 2008). R eco components respond differently to warming. Aboveground respiration component increased due to 20% increased aboveground productivity linked to soil warming of 2.3 0 C at a cold region site in Alaska (Natali et al., 2012). Root respiration showed almost no change after a 2 0 C of soil warming but a 21% increase in heterotrophic respiration was observed due to enhanced microbial activities (Wang et al., 2014). The soil carbon cycling involves multiple microbial species (Crowther et al., 2019), which require favorable redox environment in soil for their sustenance (de Angelis et al., 2010). Soil carbon mineralization has been impacted by hydrological processes too (Anthony et al., 2018). Atmospheric CO 2 concentration is reported to be influenced by terrestrial water storage change (Humphrey et al., 2018). All of the processes mentioned above are not undergoing alone, rather they are occurring with feedback from each other. Therefore, integration of hydrological, redox as well as microbial processes are necessary to enhance the model estimates (Bhanja et al., 2019a(Bhanja et al., , 2019b).
Contemporary earth system models (ESMs) employ empirical approaches (mostly as functions of either soil temperature or moisture contents) including transformation of numerous soil carbon pools at different rates for estimating R eco (Clark et al., 2011;Oleson et al., 2010;Davidson and Janssens, 2006). Widely used regional-scale models, CENTURY and DAYCENT models are using simplified functions of soil temperature and moisture to estimate soil respiration (Del Grosso et al., 2005). Similar first order kinetics-based approaches are modeled in other well known regional-scale models such as ORCHIDEE (Qiu et al., 2018), Biome-BGC (White et al., 2000) and CLM (Lawrence et al., 2019). However, the relationship between incident soil temperature and soil carbon transformation is not straightforward because microbes play crucial roles in litter and soil organic matter transformation (Melillo et al., 2017). As a result, most ESMs are unable to evaluate the effect of other environmental factors, such as dissolved oxygen and redox potential, for realistic, grid-scale estimates of soil carbon (Todd-Brown et al., 2013). Large amount of errors are found in global-model generated R s in arctic and subarctic regions with low R s values (Hursh et al., 2017). Challenges to develop such regional-scale, integrated models are twofold. Firstly, the model must be well evaluated and operated across a range of domains (e.g. climate, soil, land use, topography, hydrology) as well as land use change and management. Particularly, in cold regions, permafrost and freeze-thaw cycles due to climate change may cause a substantial change in hydrology and microbial activities, leading to change of soil respirations and soil gas exchange (Cui and Wang, 2019). Secondly, to calibrate and drive the models, various datasets with fine spatial resolution and continuous measurements at a range of different spatial scales are required. Particularly, at a regional scale, it is pertinent that the data requires to cover the diverse meteorological, soil, land use and management domains, all with the quantified uncertainties.
Keeping these issues in mind, here we introduced detailed litter decomposition, root respiration and aboveground respiration modules in order to simulate realistic estimates of R h , R a and R eco . In the recent years, Bhanja and Wang (2020) introduced the soil CO 2 emission capabilities in the Soil and Water Assessment Tool (SWAT). Indeed, Bhanja et al. (2019a) and Bhanja et al. (2019b) uniquely introduced multiple major chemical reactions and SOM decomposition, which integrated main environmental factors with the heterogeneity of vegetation, soil and hydrological processes. This enables us to evaluate the effect of environmental factors on R h , R a and R eco . Here the influence of environmental drivers on R a , R s and R eco are investigated for the first time in boreal forest covered Athabasca river basin (ARB), Canada. We perform a dominance analysis using a multivariate approach to find out the influence of individual environmental factors on R a , R s and R eco . The contribution of the 6 predictor variables (air temperature, precipitation, soil temperature, WFPS, dissolved oxygen level and ORP) on respiration rates are estimated based on the R 2 values originated from multiple regression analyses. We also attempted to validate the respiration estimates on comparing with the site-scale measurements available.

Modeling framework
The SWAT has been widely used for regional-scale simulation capacity of detailed hydrology (Arnold et al., 1998;Neitsch et al., 2011). We added additional carbon cycle capabilities into SWAT for simulating respiration components originating from litter decomposition, root, above-ground biomass and soil organic carbon transformation (Fig. 1). The new integrated hydro-biogeochemical model (named as SWAT-MKT) is capable of performing variety of tasks at one run (for detailed descriptions of the base model, please refer to Bhanja et al., 2019aBhanja et al., , 2019bBhanja and Wang, 2020). Major chemical processes are considered in this approach, as shown in Figure S1.
Different remote sensing as well as modeled datasets are used to set up the SWAT model simulation for the ARB. km × 1 km) was accessed from the United States Geological Survey (USGS; Loveland et al., 2000). Soil map (1:1 million resolution) was obtained from the Agriculture and Agri-Food Canada (AAFC; SLC, 2010). Watershed delineation was done through SWAT at a 200 km 2 threshold that characterizes 131 sub-basins in ARB ( Figure S2). The functional unit in SWAT were characterized as the Hydrologic Response Units (HRU). Four slope classes (5%, 10%, 15% and 20%) with 10%, 5% and 10% thresholds were used for land-use, soil and slope, respectively. Finally, 1370 HRUs were characterized in ARB. Daily-scale precipitation, minimum and maximum temperature data were obtained from the GoC (2016) database. Daily-scale wind speed, relative humidity and solar radiation data were accessed for 230 stations from CFSR (2016) archive. Further details on model built-up were provided in Shrestha et al. (2017).

Litter decomposition
Liter decomposition sub-module was developed assuming the existence of two soil carbon pools, active and passive, respectively. Active litter fraction along with the microbial biomass were considered in the active pool (Fujita et al., 2014). We used CENTURY model's first order litter decomposition kinetics approach as a baseline (Parton et al., 1987(Parton et al., , 1994(Parton et al., , 2001Fujita et al., 2014): Where, R di,C (gC kg − 1 soil d -1 ) is the litter decomposition rate from the CENTURY model, C i the carbon contents within active or passive substrate (gC kg − 1 soil), k i,C the first-order decomposition coefficient of Ci (d -1 ), and i the active and passive substrate type.
As microbes are playing major role in litter decomposition, microbial enzymatic approach adopting Michaelis-Menten kinetics (Fujita et al., 2014) was used here. The modified decomposition rate R di,M (gC kg − 1 soil d -1 ): Where, the half-saturation constant or Michaelis-Menten constant is represented as Km i (gC kg − 1 soil). The decomposition coefficient of C i is k i,M and it is estimated individually corresponding to the active (AC) and passive (PA) substrates: Where, decomposition coefficients used in CENTURY model for the active and passive substrates are represented ask AC,C and k PA,C , respectively. Km AC is the Michaelis-Menten constant for active substrate and it is estimated at 0.3 g C kg − 1 soil (Allison et al., 2010). Km PA is the Michaelis-Menten constant for the passive substrate and its value is estimated at 600 g C kg − 1 soil (Allison et al., 2010). Microbial biomass carbon is represented as C b (g C kg − 1 soil), its value was estimated as the global median microbial biomass (0.87 g C kg − 1 soil; Cleveland and Liptzin, 2007). Total carbon stock of soil, C T is considered as the global total soil carbon: 46 g C kg − 1 soil (Cleveland and Liptzin, 2007).
While the microbes facilitating the litter decomposition rates, its own biomass can be well subjected to active decomposition. Therefore, microbial biomass is considered as an active litter component and its magnitude is directly proportional to the litter decomposition rate. The new rate (R di,MM ) can be estimated as (Fujita et al., 2014): Soil respiration rates (R LD ) associated with litter decomposition was estimated following Fujita et al. (2014).
Where, the growth efficiency of microbes during decomposition of either active or passive substrates is represented as e i,m (0.45; Fujita et al., 2014). Overflow of carbon due to low nitrogen concentration is represented as O m,c ; we have omitted this parameter due to data limitation to represent the processes. Finally the Eq. (6) becomes:

Root respiration
SWAT model is not capable of simulating the root respiration (R r ). We have introduced a new sub-module in SWAT to simulate R r following Li et al. (1994): Where, CO 2 produced by roots for nitrogen uptake is represented as R n (13.8 mg C meq -1 N; Veen, 1981;Li et al., 1994). U n (kg N ha − 1 d -1 ) represents the nitrogen uptake rates of plant. CO 2 produced as a function of root growth is represented as R rg (19.19 mg C g − 1 dry matter; Veen, 1981;Li et al., 1994). Root biomass growth per day is represented as BG r (g dry matter ha − 1 ). CO 2 produced due to root maintenance is represented as R rm (0.288 mg C g − 1 dry matter d -1 ; Veen, 1981;Li et al., 1994). B lr represents the living root biomass (g dry matter ha − 1 ).

Respiration from aboveground biomass
Respiration by above ground biomass (R abv ) was not present in SWAT model. It is estimated following Ryan et al. (1994): Where, daily aboveground foliar biomass growth: BG abvf (g dry matter ha − 1 d -1 ). CO 2 produced due to aboveground foliar biomass growth: R abvf (1.767 mg C g − 1 dry matter d -1 ; Ryan et al., 1994). Daily aboveground woody biomass growth is represented as BG abvw (g dry matter ha − 1 d -1 ). CO 2 produced due to aboveground woody biomass growth is represented as R abvw (0.12 mg C g − 1 dry matter d -1 ; Ryan et al., 1994).

Dominance analysis
Dominance analysis was performed to find a qualitative relation defined in a pairwise fashion (Budescu, 1993). If one variable is more useful than its competitor in all subset regressions, it is called to dominate another. We adopted a multivariate approach for dominance analysis to rank order of individual environmental factors that influence on R a , R s and R eco in terms of their relative importance. The contribution of the 6 predictor variables (air temperature, precipitation, soil temperature, WFPS, dissolved oxygen level and ORP) to the respiration are estimated based on the statistical performance of R 2 values originated from multiple regression analyses (Azen and Budescu, 2006).

Assumptions and limitations
Geological CO 2 emission from mineralization was not considered in this study (Andrews and Schlesinger, 2001). Animal respiration is not included in ecosystem respiration computation. However, in this prevailing cold climatic conditions, number of animals residing at the ARB is very low (Weber et al., 2015), respiration from animals can be neglected here. Several assumptions and limitations related to the basic version of the model were described in Bhanja et al. (2019a) and Bhanja et al. (2019b) and Bhanja and Wang (2020). All of the soil microbes are either not participating in the decomposition process or even participating and performing decomposition at a slower rates leading to restriction of the decomposition process (Kaiser et al., 2015). Values of some of the constants in this respiration modeling approach are taken from literature information due to non availability of the data within the study region.

Trends of environmental factors, total ecosystem respiration and its components
Annual mean air temperature, soil temperature and precipitation show spatial variations across ARB ( Figure S2, S3, S4). Annual trends of air and soil temperature show an increasing pattern during 2000-2013 in most of the subbasins of ARB (Fig. 2a, c). Precipitation trends show a distinct north-south demarcation with increasing values in southern subbasins and decreasing values in northern subbasins (Fig. 2b). WFPS trends show a mixture of increasing and decreasing patterns while dissolved oxygen levels show an increasing trend in most of the subbasins (Fig. 2d, e). WFPS as a measure of soil water availability, it is related to precipitation, temperature, snow melt as well as the freeze and thaw cycle of the soil (Bhanja and Wang, 2021). WFPS is also dependent on the soil physical properties that are not homogenous in the study region (Bhanja et al., 2019a). These are the reasons for its distinct deviation on comparing with precipitation patterns. Dissolved oxygen level is dependent on multiple factors, such as incident oxygen concentration, nature of chemical reactions, air filled pore space and soil temperature (Bethke, 2007;Fan et al., 2014). ORP trends show increasing values in all of the subbasins with a couple of exceptions (Fig. 2f). Soil ORP being the function of oxygen concentration, microbial activities and the prevailing chemical reactions (Bethke, 2007;Reddy and Delaune, 2008;Bhanja and Wang, 2020), it is very difficult to state the exact reason for its rise over the years. The ORP increase would make the system more oxidized and hence if it continues to rise in future, soil chemical dynamics would change. Earlier studies used WFPS as a proxy of ORP to estimate soil greenhouse gas emissions (Bateman and Baggs, 2005;Wagena et al., 2017;Yang et al., 2017;Shrestha et al., 2018), however, differences in WFPS and ORP trends are distinct. R a , R s and R eco trends show near similar spatial patterns of increasing and decreasing values but the magnitudes are different (Fig. 2g, h, i). As observed for trends of temperature (both air and soil) and precipitation, temperature shows rising trend and precipitation shows declining trend in most of the northern subbasins. The spatial patterns of temperature and precipitation trends are not evident for R a , R s and R eco in northern subbasins. The modeled estimates compare well with the available R eco estimates (Table S1) at nearby boreal sites taken from FLUXNET measurements (Pastorello et al., 2017). Simulated R s values are also matching well with nearby site data (Table S2) available from SRDB archives (Bond-Lamberty and Thomson, 2010b). It should be noted that the majority (if not all) of the FLUXNET and SRDB sites are located in the south of the ARBthis is the reason for the observing some higher values at the sites.

Influential factors of autotrophic and soil respirations
We have investigated the relationship between R a and the environmental influencing factors (Fig. 3). Temperature (both air and soil) variations can explain R a variations well at most of the subbasins as statistically significant (p < 0.001) correlations are observed at most of the subbasins of ARB (Fig. 3a, 3c). Dissolved oxygen and ORP values can also explain R a values well with very good correlation coefficients (Fig. 3e, 3f). Correlation of R a with precipitation and WFPS are not so good even mostly negative r values are obtained for WFPS (Fig. 3b, d). Based on the correlation analyses, soil temperature, dissolved oxygen level and air temperature are the three most crucial parameters that can well explain the R a variations. Monthly basin averaged values of R a show non-linear relationship with environmental drivers ( Figure S5). Coefficients of determination of these non-linear relationships show best results for ORP, soil temperature and dissolved oxygen level, respectively. Dominance analysis for the monthly-scale basin-averaged values indicated complete dominance of soil temperature on influencing R a with 24% contribution (Fig. 4). Dissolved oxygen and air temperature are the next two most influencing factors of R a .
Prevailing temperature (air and soil both) influence the R s most with observed correlation coefficients > 0.9 in most of the subbasins for soil temperature (Fig. 5a, c). Dissolved oxygen level and precipitation are the next two main influential drivers of R s (Fig. 5b, e). ORP closely follows precipitation in most of the subbasins or even performs better in    a handful of them (Fig. 5b, f). WFPS found to be the least influential factor for R s (Fig. 5d). Scatter analysis using basin-averaged environmental factors exhibits non-linear relationships with R s ( Figure S6). The coefficient of determination is found to be best for soil temperature, air temperature and dissolved oxygen level. Dominance analysis also shows strong influence of soil temperature, followed by dissolved oxygen and air temperature on R s (Fig. 4). The three parameters contribute nearly 67% of R s . Influence of temperature on R s is well reported in earlier studies (Bond-Lamberty and Thomson, 2010ab;Chen et al., 2010;Hursh et al., 2017). Our observation based on dissolved oxygen is not reported before.

Influential factors of ecosystem respiration
Soil temperature and dissolved oxygen levels are influencing R eco most with statistically significant (p < 0.001), very high correlation (r > 0.8) at most of the subbasins (Fig. 6c, e). Air temperature, ORP and precipitation are the other important drivers of R eco with moderate to good correlation estimates (Fig. 6a, b, f). Similar to other respiration estimates, WFPS and R eco relationships are weak (Fig. 6d).
Seasonal cycles are distinctly observed in basin-averaged R eco and other environmental drivers ( Figure S7). ARB has been characterized by cold and dry winter and wet summer ( Figure S7), which is the reason for the seasonal cycle. Our model performs well on comparing the simulated WFPS and soil temperature with the observed values at 3 locations (Bhanja et al., 2019a). Temperatures (both air and soil), precipitation and dissolved oxygen level show lowest values during winter time. WFPS exhibit maximum values during early summer/snow melt period and lowest values during late summer period when evapotranspiration rates are higher. ORP shows lowest values in late winter during March-April ( Figure S7).
Scatter analyses of basin averaged R eco and environmental factors show nonlinear patterns (Fig. 7). Soil temperature, dissolved oxygen levels and air temperature are ranked first three variables based on coefficient of determination of these non-linear fitting patterns. Based on the scatter analyses, WFPS is not capturing the R eco patterns well (Fig. 7). Dominance analysis based estimates strongly support our findings from correlation and scatter analyses (Fig. 4). Out of the 6 environmental drivers we have used, soil temperature completely dominated R eco with nearly 26% contribution while the dissolved oxygen and air temperature are the next two influencers of R eco with estimated contributions of 20% and 19%, respectively (Fig. 4).

Conclusions
We have studied the influence of environmental factors on R a , R s and R eco using our newly developed model at the Athabasca River basin, Canada. Our model with detailed microbially-mediated, carbon cycle capability enables users to simulate R a , R s and R eco at a daily time step. A dominance analysis was performed using a multivariate approach to rank the influence of the 6 predictor variables (air temperature, precipitation, soil temperature, WFPS, dissolved oxygen level and ORP) on R a , R s and R eco, based on the R 2 values originated from multiple regression analyses. Respiration estimates show clear spatial as well as temporal variations. Our results show that the air and soil temperature and dissolved oxygen level are exhibiting the highest correlations at most of the subbasins on comparing with all of the respiration estimates and are the major influencing factors with >65% contribution for the respiration estimates (R a , R s and R eco ). Scatter and time-series analyses also support the result.
In general, the previous studies reported either temperature and soil moisture or both in combination as the major influencing factors of R eco and R s . Results from this study indicate other environmental drivers such as dissolved oxygen levels and the ORP of the system play a crucial role in R eco dynamics. More studies using site-scale data remain to be required to improve our understanding in this rarely studied area. We believe incorporation of these components in future models could improve regional-scale respiration simulation.

Credit statement
S.N. Bhanja conceived the idea of the manuscript, collected the data, developed the model and conducted the statistical analyses. S.N. Bhanja has written the manuscript with inputs from J. Wang.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.