Evaluating an exponential respiration model to alternative models for soil respiration components in a Canadian wildfire chronosequence (FireResp, v1.0)

Forest fires modify soil organic carbon and suppress soil respiration for many decades since the initial disturbance. The associated changes in soil autotrophic and heterotrophic respiration from the time of the forest fire however, is less well characterized. The FireResp model predicts soil autotrophic and heterotrophic respiration parameterized with a novel dataset across a fire chronosequence in the Yukon and Northwest Territories of Canada. The dataset consisted of soil incubation experiments and field measurements of soil respiration and soil carbon stocks. The FireResp model contains submodels that 5 consider a Q10 (exponential) model of respiration to models of heterotrophic respiration using Michaelis-Menten kinetics, parameterized with soil microbial carbon. For model evaluation we applied the Akaike Information Criterion and compared predicted patterns in components of soil respiration across the chronosequence. Parameters estimated with data from the 5 cm soil depth had better model-data comparisons than parameters estimated with data from the 10 cm soil depth. The model-data fit was improved by including parameters estimated from soil incubation experiments. Models that incorporated microbial carbon 10 with Michaelis-Menten kinetics reproduced patterns in autotrophic and heterotrophic soil respiration components across the chronosequence. Autotrophic respiration was associated with aboveground tree biomass at more recently burned sites, but this association was less robust at older sites in the chronosequence. Our results provide support for more structured soil respiration models than standard Q10 exponential models.

In high-latitude forests, soil respiration fluxes and soil carbon stocks exhibit variation depending on the time since the last 30 wildfire (Bond-Lamberty et al., 2004;O'Donnell et al., 2011). Fire modifies soil organic carbon quality, making it harder for microbes to access carbon (Holden et al., 2016;Song et al., 2019;Zhao et al., 2021). A recent meta-analysis by Ribeiro-Kumara et al. (2020b) of 32 studies measuring soil respiration following wildfires indicates two emergent patterns. First, overall soil respiration stabilizes 10-30 years following a fire. Second, for components of soil respiration, R A will increase and ultimately approach a steady-state value associated with forest succession and vegetation regrowth. On the other hand R H may decrease 35 by association with post-fire changes in soil organic matter quality, temperature, or moisture (Aaltonen et al., 2019a, b;Wei et al., 2010). For a sense of the magnitude of these changes, Bond-Lamberty et al. (2004) found the proportion of annual soil respiration that is R A changes from 5% (following disturbance), to 40% (21 years post-disturbance), and returning to 15% (150 years post-disturbance). The robustness of any patterns in R A and R H is highly uncertain given known soil heterogeneity in these high-latitude soils (e.g. permafrost versus non-permafrost soils, microbial versus fungal species composition).
2. When corroborated with ::::: tested :::::: against : observational data, soil models that incorporate microbial carbon will better replicate the observed dynamics and associated fluxes (R A , R H , and the ratio R A /R S ) across the fire chronosequence.
To evaluate our hypotheses we combine data from soil incubation experiments (Aaltonen et al., 2019b) with field data (Köster et al., 2017) at chronosequence sites. For both incubation and field data, measurements were collected at the same time from similar plots to minimize any spatial and temporal biases in the data. Models :::::::: Submodels : are evaluated based on their ability to 70 replicate measured soil respiration (both from incubation and field measurements). To reduce any biases with model fitting or model equifinality (Christiansen, 2018;Marschmann et al., 2019) we evaluate a range of parameter estimation approaches and data types.  et al. (2017) and Aaltonen et al. (2019b). 80 Chronosequence sites were selected from the time since last burned with a stand replacing ::: the ::: last : fire (in 1968, 1990, and 2012) and ::: that :::::: burned :: all ::::::::::: aboveground ::::::::: vegetation. ::: We :::: also included a control site, where the last fire was more than 100 years ago. The date and boundaries of the fires were determined from geographic data from the Canadian Wildfire Information System (Natural Resources Canada). We visually corroborated the geographic location of our sites with reported fire boundaries.

110
Incubation data also measured :::: The ::::::::: incubation :::: data ::::::: included :::::::::::: measurements ::: of the available soil organic carbon extracted from incubation soils, denoted here as C A , as described in Zhou et al. (2019). Briefly, soil dissolved organic C content was measured using total organic C analyzer (Shimadzu TOC-V CPH, Shimadzu Corp., Kyoto, Japan) from soil extracts extracted with 0.5 M K 2 SO 4 . Microbial carbon used in the FireResp model was extracted using the chloroform fumigation extraction method (Beck et al., 1997). Briefly, three grams dry weight equivalent of soil was fumigated at 25 • C with ethanol-free 115 chloroform for 24 h and extracted with 0.5-M K 2 SO 4 . The conversion factors, also known as the extraction efficiency, for estimating the microbial carbon is 0.45 (Beck et al., 1997). For the field data, we approximated C A as linearly associated with total soil carbon C S at a given depth, extrapolated from linear regression in the incubation data (results not shown).
For the field samples an estimate of root carbon C R was assumed to be proportional to total tree biomass collected at each plot (Härkönen et al., 2011;Neumann et al., 2020). A summary of all input variables is reported in Table 1.  Table 1. Summary of soil measurements for this study, organized by site in the chronosequence ( , 1990( , 1968 and depth of measurement. ::: The :::
(1) temperature is a commonly applied (empirical) paradigm for respiration, motivated by temperature dependencies of enzymatic reactions (van't Hoff and Lehfeldt, 1898). This exponential temperature model is applied for R A and R M , similar to process models for these components at the ecosystem scale (Aber et al., 1997;Zobitz et al., 2008). Input variables for 130 Eq.
(2) arises from first-order microbial enzyme kinetics (Allison et al., 2010) under quasi-steady state assumptions (Keener et al., 2009). In Eq. (2), is the efficiency converting substrate to microbe :::::::: microbial : biomass (no units), µ is the maximum microbial uptake rate (hr −1 ), and k A (g C m −2 ) represents the half-saturation rate, and C X represents the substrate for respira-145 tion. Depending on the model variant C X may be total soil carbon (C S ) or available soil organic carbon (C A ), which represents more labile carbon for ingestion by microbes.
(2). With these considerations total soil respiration is expressed in Eq. (4): The Microbe submodel is based on a two-pool soil-microbe model described in Sihi et al. (2016).
-Microbe-mult submodel: This submodel is structured similar to the Microbe model but with two modifications. First, growth respiration is not considered. Second, maintenance respiration is multiplied by a Michaelis-Menten factor: The Microbe-mult model is designed to be an intermediate model between the Null model and the Microbe model. The

165
additional multiplicative factor is a heuristic designed to represent maintenance respiration as substrate limited by C S .
-Quality submodel: This submodel is structured similar to the Microbe model, but for growth respiration (R G ) available soil organic carbon (C A ) is the input for pool C X in Eq.
(2). Total soil respiration is expressed in Eq. (6): The Quality submodel is based on a multi-pool soil model that structures the soil into different pools based on the 170 recalcitrance and turnover time of the soil parent material, similar to models by Bosatta and Ågren (1985). Inputs from litterfall, enzymatic degradation, root turnover, or root exudation create a pool of available soil organic carbon (C A ) that can be incorporated into microbial biomass. While in this case R G is represented with Eq.
(2), when constructing a dynamic model of soil would additionally include expressions for the transformation of each soil pool through enzymatic degradation and mineralization to a more recalcitrant pool (both under first-order kinetics).

175
-Quality-mult submodel: This submodel is structured similar to the Quality model with two modifications (similar to the modifications made in the Microbe-mult model). First, growth respiration is not considered. Second, maintenance respiration is multiplied by a Michaelis-Menten factor: Like the Microbe-mult model, the Quality-mult is a heuristic model designed to represent maintenance respiration as 180 substrate limited by C A . the ::::::: FireResp : model along with the allowed range. * : denotes a parameter for the Incubation Field Linear parameter estimation approach. †: denotes a parameter for the Field Linear parameter estimation approach. Table 2 summarizes the different parameters for each model and their allowed ranges when estimating parameters.

Parameter estimation routine
The different submodels (Null, Microbe, Quality, Microbe-mult, and Quality-mult) may be nonlinear with respect to the parameters. For parameter estimation we applied the Levenberg-Marquardt algorithm (Elzhov et al., 2016). The Levenberg-Marquardt 185 algorithm optimizes an objective function, which in this case is the residual sum of squares between measured and modeled soil respiration R S . The algorithm also requires (1) the Jacobian of the model to accelerate convergence to the optimum value, (2) an initial guess for parameters, (3) and bounds for all parameters.
The Levenberg-Marquardt algorithm may converge to a local (rather than global) optimum or the estimated parameter values may be at the boundaries of the allowed range. To ensure that parameter estimates converged to a global (rather than 190 local optimum), initial parameter guesses for the method were drawn from a uniform distribution with reasonable bounds on parameters ( Table 2). The Levenberg-Marquardt algorithm is implemented in R with the package nlsr (Nash, 2014;Nash and Murdoch, 2019).
(1) Field: All model parameters (e.g. Q 10,M , k M , k A , µ, , k S , Q 10,R , and k R ; depending on the type of model) were estimated with the field data only.
(2) Field Linear: Model parameters for R H (e.g. Q 10,M , k M , k A , µ, , k S , depending on the type of model) are estimated with the field data. Rather than a Q 10 function for R A (Eq. (1)), for this approach R A equals g R · C R , where C R is provided by the field data. We then estimated g R from the field data.

205
(3) Incubation Field: Two separate parameter estimations were applied. First, : model parameters for R H (e.g. Q 10,M , k M , k A , µ, , k S depending on the type of model) were first estimated with the incubation data. Next, autotrophic respiration parameters (Q 10,R and k R ) were estimated from field data.
(4) Incubation Field Linear: Similar to the Incubation Field approach, parameters relating to R H were first estimated with incubation data. Next using these parameter estimates, heterotrophic respiration was computed from the corresponding 210 field measurements (denoted here as R H,f ield ). Total soil respiration then equals R S = g R · C R + f · R H,f ield , with R A = g R · C R and R H = f · R H,f ield . We then estimated f and g R from the field data. Table 3 shows the relationship between the different parameter estimation approaches studied. Table 4 lists the parameters estimated for each model :::::::: submodel : and parameter estimation approach. Data used for parameter estimation consisted of combinations from five different categories of sites ( , 1990( , 1968 or all sites 215 together) and 3 different depths (5 cm, 10 cm, or both depths together). Additionally with the four different parameter estimation approaches (Field, Field Linear, Incubation Field, and Incubation Field Linear) and five different models ::::::::: submodels (Null, Microbe, Microbe-mult, Quality, and Quality-mult), 300 separate parameter estimations were computed.

Model evaluation
We applied two different approaches to evaluate the reasonableness of a model-data fit. The first approach relied on Taylor diagrams (Taylor, 2001), which facilitates intercomparison between models when compared to measured values (in this case R S ). The Taylor diagram is structured as a polar coordinate plot; here the radius ν is the normalized ratio between modeled 230 and measured standard deviation σ model /σ measured and the angle θ corresponding to the correlation coefficient r for measured and modeled R S . Two comparisons can be visually inferred from the Taylor diagram. First, the point located at (ν, θ) = (1, 0) represents a set of modeled values of R S that perfectly match measured R S . Values of ν less than unity indicate that modeled R S has less variability. Second, the distance from a point on the diagram to (ν, θ) = (1, 0) is the centered pattern root mean square distance. Concentric circles from the point (ν, θ) = (1, 0) help assess the compare the centered pattern root mean square 235 distance for modeled results.
A second approach relies on Akaike's Information Criteria (AIC) (Akaike, 1974). The AIC is defined as −2·LL+2·p, where LL is the log-likelihood and p the number of parameters in the model. The submodel with the lowest AIC is defined as the best approximating model for the data. We apply the AIC to compare across submodels for a parameter estimation approach to control for sample size effects in the AIC.
Columns in the facetted plot represent the depth of the data used for parameterization (5 cm, 10 cm, or All depths), rows refer to the chronosequence site. Radii represent the normalized standard deviation between a :::::: FireResp : submodel value of RS to measured RS, angles the correlation coefficient r (labeled). The dashed concentric circles represent contours (increments 0.25) for the normalized centered pattern root mean square distance.
the range of the parameters :::::: allowed ::::::::: parameter ::::: range) as separate colors. In contrast, Figure 6 structures each sparkline plot by the submodel studied (Null, Microbe, Quality, Microbe-mult, and Quality-mult), facilitating comparisons between models for a given parameter estimation and depth of data used in the parameter estimation. In Figure 6, sparkline plots for adjusted R 2 or 260 AIC values are all scaled respectively the same for each statistic. The models with the largest adjusted R 2 or lowest AIC value are denoted as separate colors.
We computed R A , R H , and the proportion of soil respiration due to autotrophic respiration (p A = R A /(R A +R H )) for each parameter set generated through the parameter estimation routine (Section 2.3). We then computed summary statistics from the distribution of R A , R H , p A for each parameter estimation approach. Summary results for the median of these distributions 265 for R A and R H are shown in Figure 7, organized by the parameter estimation approach. Additionally the red shading in Figure 7 shows the minimum and maximum ranges of measured R S (lines), first or third quarters (boxes), and median R S for comparison. Figure 7 visually displays no significant difference in patterns of R A and R H by the depth of the soil data used for parameter estimation (5 cm, 10 cm, or both depths together). (2020a)). We computed the predicted values of p A from a loess fit using years since disturbance and p A as variables. Figure 5. Median values of parameter estimates for different models ::::::: FireResp :::::::: submodels using the Incubation Field Linear approach at 5 cm cm depth. The horizontal axis on each sparkline plot is arranged by the year since the burn sites in the chronosequence ( , 1990( , 1968. In each column the vertical axis scale is the same. Edge hitting ::::::::: Edge-hitting : parameters :::::: (defined ::: here :: as :::::: within : a :::: tenth :: of :::::: percent : of ::: the :::::: allowed :::::::: parameter ::::: range) are denoted with the blue coloring. Figure 6. Median values of the adjusted R 2 and AIC from different parameter estimation approaches (Field, Field Linear, Incubation Field, and Incubation Field Linear) using measurements made at given depth. The horizontal axis on each sparkline plot is arranged by the different models studied a :::::::: FireResp :::::::: submodels (Null, Microbe, Quality, Microbe-mult, and Quality-mult). For the adjusted R 2 sparkline plot, the vertical axis ranges between 0 to 1, with gridlines every 0.25 units. The submodel with the highest adjusted R 2 value is denoted with red coloring. For the AIC plots, the vertical axis ranges from 50 to 150 ::: 250, with gridlines every 50 units. The submodel with the lowest AIC is denoted with red coloring. Incubation Field Linear), soil depth data used for parameter optimization (5 cm, 10 cm, or both depths together) and submodel (Null, Microbe, Quality, Microbe-mult, and Quality-mult). The ::: grey :::: lines ::: are :::: used :: as : a ::::: guide :: to :::: show ::: the :::::::::::: chronosequence :::: trend ::: for : a :::::::: particular ::::::: parameter :::::::: estimation ::::::: approach ::: and ::: soil ::::: depth. ::: The : boxplot shows measured ranges of RS at each site in the chronosequence. Soil models that directly incorporated microbe ::::::: microbial : carbon produced patterns of R A and R H that increased from the time 275 since the fire ( Figure 7). As these patterns also conform to changes in root carbon (which was proportional to tree biomass, Table   1), we have initial support for our two primary hypotheses: (1) autotrophic respiration should be positively associated with the time since disturbance because of changes in aboveground foliar vegetation from forest succession and (2) when corroborated with ::::: tested :::::: against : observational data, soil models that include soil microbial carbon will better replicate expected patterns for soil respiration components across the chronosequence. We will further evaluate the two hypotheses through subsequent 280 analysis of the data used for parameter estimation, parameter estimation approaches, and the soil respiration models.

Evaluation of datasets for parameter estimation
We had two categories of datasets for this study: the type of data (incubation or field data) or the depth at which measurements were made (5 cm, 10 cm, or both depths together). This controlled experimental design is also represented in the Taylor diagrams ( Figure 3) which shows, comparatively, a centered pattern root mean square distance (distance between a point on 285 the Taylor diagram and (ν, θ) = (1, 0)) ranging from 0.25 -1 and r ranging 0.7 -0.9. For the field data (Figure 4), the centered pattern root mean square distance ranged from 0.5 -1 and r 0.4 :: 0.3 : -0.8 ::: 0.9. We attributed the difference ::::::::: differences between Figures 3 and Figure 4 is that the range of soil temperatures from the incubation experiments spanned 1 -19 • C, allowing for a wider :::::::::: temperature range to characterize any exponential temperature profile. In contrast, field measurements ranged from 4 -9 • C (Table 1). For both Figure 3 and Figure 4, : the 5 cm depth had higher values for r and a smaller centered pattern root mean 290 square distance compared to the 10 cm depth.
We did not find any noticeable site differences in submodel outputs depending on the depth of the soil used for data assimilation (5 cm, 10 cm, or both depths together; Figures 3, 4, Figure 6). While soil model parameters (such as Q 10 ) are expected to vary with soil depth (Pavelka et al., 2007;Graf et al., 2008;Pumpanen et al., 2008) we did not observe any significant depth-dependent differences in parameter estimates (see the figures in the Supplementary Information). The primary reason 295 for this result is that the inter-site variability is larger than the variability by depth at a given site (Table 1 and Figure 2). We also did not find any improvements in our results when all data from sites were pooled together ( Figure 7 and Figure 8). From these conclusions : , we will limit the discussion to evaluating model results generated from data at the 5 cm depth.

305
While there is no universal pattern to R H following forest fire disturbances (Ribeiro-Kumara et al., 2020b), we have reason to believe the near-zero modeled values for R H for the 1968 site in Figure 7 may be an underestimate. For our sites we expect modest, and perhaps decreasing (but not zero), changes in R H from the time of disturbance due to three reasons. First, factors influencing recovery of R H are burn severity or intensity (Meigs et al., 2009;Hu et al., 2017) and decomposition of pyrogenic litter (Kulmala et al., 2014;Muñoz-Rojas et al., 2016). The fires at our sites combusted a significant amount of soil 310 organic matter (Köster et al., 2017) resistant to decomposition (Knicker, 2007;Aaltonen et al., 2019a), thereby minimizing any increases in R H from the decomposition of labile litter. Additionally, from this chronosequence, Aaltonen et al. (2019b) reported increased temperature sensitivity (Q 10,M ) in recently burned sites, but this was tempered by decreases in soil organic matter quality (Aaltonen et al., 2019a). Second, as succession occurs, the increase in aboveground vegetation insulates the soil, decreasing the active layer and thereby decreasing R H (Köster et al., 2017). Third, at the same chronosequence sites Our models implicitly assumed an increasing exponential relationship between temperature and respiration. The temperature sensitivity of respiration (Q 10 ) across ecosystems can vary (usually around 2-5) (Chen and Tian, 2005;Wang et al., 2006;320 Bond-Lamberty and Thomson, 2010) and is generally expected to be greater than 1, but the Q 10 value may decrease as soils warm (Niu et al., 2021). Some degree of additional variability is expected when considering the biochemical or thermodynamic foundations of respiration (Lloyd and Taylor, 1994;Ito et al., 2015), the methodological approach used to measure soil respiration (Ribeiro-Kumara et al., 2020b), or variation in the soil organic matter supply (Davidson et al., 2006).
However, an increasing exponential relationship between temperature and respiration may not be robustly supported with 325 observed data at the chronosequence sites. The forest fires at each site burnt a large portion of soil organic matter and killed the roots. Immediately following a fire, R S will be lower even if there are higher soil temperatures. In late-successional forests, the soil is colder and the active layer depth is smaller, even though there may be more soil respiration due to higher amounts of roots and soil organic matter; we observed such patterns across the chronosequence. The 2012 and 1990 sites had the highest values of T soil (Table 1) but the lowest overall respiration (Figure 7). Across the chronosequence, scatterplots of respiration with 330 temperature had a null or a negative relationship with temperature (results not shown). Empirically the negative association of respiration with temperature would imply a Q 10 value less than unity. As a result, to compensate for these opposing tendencies the R H parameters tend to be edge-hitting ( Figure 5 and Supplementary Information).
We recommend either the Incubation Field or Incubation Field Linear parameter estimation approach for two reasons. First, values of the proportion of the respiration that is autotrophic (p A = R A /(R A + R H ), Figure 8) for the Field or Field Linear 335 approaches are unexpectedly and unrealistically large, attributed to the variation in R H (Figure 7). As a baseline, Hanson et al. (2000) reported values of R A /(R A + R H ) to be approximately 0.50, which has also been supported in meta-analyses (Soil Respiration Database, Bond-Lamberty and Thomson (2010)). Second, the Incubation Field and Incubation Field Linear approaches in Figure 8 show a temporal pattern in p A similar to patterns reported in Bond-Lamberty et al. (2004)

Evaluation of hypotheses
Our first hypothesis concerned the dependence of R A on tree biomass. We developed this hypothesis from our previous studies, which concluded tree biomass was a key factor explaining patterns of soil respiration across the chronosequence (Köster 345 et al., 2017;Aaltonen et al., 2019a, b). For all models ::::::::: submodels : and the Field Linear or Incubation Field Linear parameter estimation approaches, R A is proportional to C R , which is proportional to tree biomass. Values of C R increase across the chronosequence (Table 1). However even with this proportional association, the results in Figure 7 indicate less support for our first hypothesis for two reasons. First, some modeled values R A at the 1990 site are higher than expected, especially given the association with R A to C R . Since C R is still comparatively low at this site, we might expect R A (and by association p A ) to be 350 near zero as well. Additionally, the near-zero values of R A are not a consequence of poorly-defined parameter estimates. The autotrophic respiration parameters are not overwhelmingly edge hitting ( ::::::::: parameters :::::: relating :: to :::: R A :::: (k R , :::::: Q 10,R , :: or :::: g R ) ::::: being :::::::: estimated :: as ::::: zero. ::::::::: (Otherwise :::: the ::::: values ::: for ::::: these :::::::::::::: aforementioned ::::::::: parameters :: in : Figure 5 or Supplementary Information), indicating appropriately defined parameter bounds in Table 2. :: the ::::::::::::: Supplementary :::::::::: Information ::: for ::::: these ::::::::: parameters ::: for ::: all ::: the ::::::: different :::::: models :::: and ::::::::: approaches :::::: would :: be ::::::::::: edge-hitting ::: and :::::::: indicated :::: with ::: the :::: blue ::::::: colored ::::: dots.) Second, and perhaps more 355 importantly, all parameter estimation approaches in Figure 7 predict R A to decrease between the 1968 and Control sites. The modeled decreases in R A are a result of observed decreases in R S (Figure 7) as C S increases. To compensatefor this : , estimated parameters k R or g R decrease across the chronosequence sites ( Figure 5 or Supplementary Information). The patterns to :: of k R or g R may be due to the parameter estimation routine compensating ::: for the confounding effects of increasing C R with decreasing R S . In summary, even though there is evidence for association between R A and tree biomass in earlier chronose-360 quence sites (2012 and 1990 sites), additional work is needed to explain reasons for the decline in R A for later chronosequence sites (1968 and Control sites). Future work could quantify field estimates of root mass, production, and turnover (Kalyn and Van Rees, 2006;Steele et al., 1997) to corroborate the values of C R used here and with the estimated decreases in k R across the chronosequence.
Our second hypothesis concerns the structural representation of soil respiration for soil models. Our submodels are arranged 365 on a continuum of complexity (Null, Microbe, Quality, Microbe-mult, or Quality-mult). When parameterizing more complex models parameters may be non-informative and/or edge-hitting (Zobitz et al., 2011). Reducing parameter dimensionality is a key consideration for model-data assimilation in the carbon cycle (Tang and Zhuang, 2008;Luo et al., 2009;Kraemer et al., 2020). Considering the Incubation Field Linear approach only, across the range of models ::::::::: submodels : the Microbe submodel had the smallest percentage of edge-hitting parameters (10%), ranging from 30 -50% for the other models.

370
While the AIC suggests a preference towards the Null submodel, we do not believe it is a sufficient criterion to choose it over the Microbe and Quality submodels. There was no noticeable improvement with the Null submodel in the Taylor diagrams for the field data (both in the values of r and the centered pattern root mean square difference. : ; Figure 4) or with the adjusted (the adjusted R 2 values in Figure 6 ranged from from .25 -.61), no submodel significantly improved the adjusted R 2 or AIC.
A design constraint was to construct models with the greatest potential to be fully parameterized from the collected data.
For the Quality-mult and Microbe-mult models ::::::::: submodels, k A was estimated at the lower end of its range ( Figure 5), essen-380 tially reducing these models to the Quality and Microbe models ::::::::: submodels respectively. Even though we cannot definitively conclude which of the two models :::::::: submodels : (Quality or Microbe) is the better approximating model, we recommend that some consideration of microbial growth and maintenance respiration be considered using Michaelis-Menten kinetics as a starting point (Davidson et al., 2006). Several frameworks already exist for incorporating Michaelis-Menten kinetics (Todd-Brown et al., 2012) or substrate quality degradation Ågren, 1991, 2002). Continuous (daily or sub-daily) soil respiration 385 measurements could better support more complex soil models (Rayment and Jarvis, 2000;Subke et al., 2006;Subke and Bahn, 2010;Phillips et al., 2011;Pumpanen et al., 2015;Zhang et al., 2015). Each of the models could be incorporated into a dynamic model of ecosystem carbon cycling (Zobitz et al., 2008) that also include temporal changes of permafrost active layer depth (Zhu et al., 2019).

390
We examined the ability to parameterize a range of soil respiration models using data collected from a fire chronosequence.
Importantly : , we found support for parameterizing a more complex submodel to replicate patterns in soil respiration and its components across a fire chronosequence. Separate analysis of soils with incubation experiments reduces the number of parameters to be estimated, however care must be taken in scaling incubation studies to field measurements.
For these high-latitude sites, future work could couple the models here to more continuous measurements of soil temperature 395 along with a dynamic active layer depth model (Zhu et al., 2019). These modeling approaches could examine the effects of gross primary productivity on soil respiration components (Zhuang et al., 2002;Pumpanen et al., 2003;Vargas et al., 2010;Pumpanen et al., 2015;Phillips et al., 2017). For sites that cannot be instrumented continuously (such as the ones studied here), this model data integration could be supported with periodic surveys of aboveground biomass and other remote sensing data (Neumann et al., 2020).

400
Code and data availability. Code and data necessary to reproduce all results are available through GitHub https://github.com/jmzobitz/ FireResp and archived in Zenodo (Zobitz et al., 2021).
Supplementary Information includes sparkline parameter estimates (similar to Figure 5) for all approaches and depths examined.
collected the field data. Co-authors Aaltonen and Zhou analyzed the incubation data. All authors contributed to evaluating the results and the writing of the manuscript.} Competing interests. The authors declare no competing interests.
Acknowledgements. Co-author Zobitz was funded by the Fulbright Finland Foundation and Saastamoinen Foundation Grant in Health and Environmental Sciences. This work was funded by the Academy of Finland (project numbers 286685, 294600, 307222, 327198, and 337550).

410
Travel funding was provided from EU InterAct (H2020 Grant Agreement No. 730938). Co-author Zobitz acknowledges B. S. Chelton for helpful discussion on this manuscript.