Biome-scale temperature sensitivity of ecosystem respiration revealed by atmospheric CO2 observations

The temperature sensitivity of ecosystem respiration regulates how the terrestrial carbon sink responds to a warming climate but has been difficult to constrain observationally beyond the plot scale. Here we use observations of atmospheric CO2 concentrations from a network of towers together with carbon flux estimates from state-of-the-art terrestrial biosphere models to characterize the temperature sensitivity of ecosystem respiration, as represented by the Arrhenius activation energy, over various North American biomes. We infer activation energies of 0.43 eV for North America and 0.38 eV to 0.53 eV for major biomes therein, which are substantially below those reported for plot-scale studies (approximately 0.65 eV). This discrepancy suggests that sparse plot-scale observations do not capture the spatial-scale dependence and biome specificity of the temperature sensitivity. We further show that adjusting the apparent temperature sensitivity in model estimates markedly improves their ability to represent observed atmospheric CO2 variability. This study provides observationally constrained estimates of the temperature sensitivity of ecosystem respiration directly at the biome scale and reveals that temperature sensitivities at this scale are lower than those based on earlier plot-scale studies. These findings call for additional work to assess the resilience of large-scale carbon sinks to warming.


1
Combination tests to assess impact of the uncertainty in gross primary productivity on the inferred temperature sensitivities of ecosystem respiration . . . . . . . . . 21 2 Potential impact of lateral emissions on the inferred temperature sensitivities . .   Fig. 1 | Optimal North American and biome-specific temperature sensitivities of ecosystem respiration for individual models and the model ensemble based on GPP combination tests. Relationships between mean adjustments to model-specific estimates of the activation energy needed to maximize consistency with observed atmospheric CO 2 variability in GPP combination tests (Δ a , vertical axis) and the original estimates of activation energy ( a , horizontal axis) for models in the MsTMIP v2 (pink), TRENDY v6 (brown), and FLUXCOM (orange) model ensembles in (a) the North American domain, (b) croplands, (c) evergreen needleleaf forests, and (d) deciduous broadleaf and mixed forests. Similar to Fig. 2 in the main text, only models for which the explanatory power of GPP estimates ( 2 GPP ) exceeds that of shortwave radiation ( 2 SW ) are included, yielding = 29 models. For each model, the mean adjustment to the activation energy (vertical axis) is determined from averaging such adjustments from 28 "mix-and-match" pairs of GPP and ecosystem respiration estimates (see Supplementary Notes 1), with vertical error bars indicating ±1 standard deviation. The gray dashed lines represent the best ODR fit between Δ a and a estimates across the models, with light gray shading indicating the 95 % prediction interval. The optimal temperature sensitivity corresponds to the point where the ODR fit line crosses Δ a = 0 eV (i.e., no adjustment to a is needed) and is indicated by a green dashed line. GPP, gross primary productivity; ODR, orthogonal distance regression.  Fig. 1, only models for which the explanatory power of GPP estimates exceeds that of shortwave radiation are included ( = 29). Vertical error bars indicating ±1 standard deviation among the adjustments derived from 28 "mix-and-match" pairs of GPP and ecosystem respiration estimates (see Supplementary Notes 1). Gray dashed lines are the best fit lines from ordinary least-squares linear regression, with light gray shading indicating the 95 % prediction interval. GPP, gross primary productivity.  Fig. 2 in the main text, only models for which the explanatory power of GPP estimates ( 2 GPP ) exceeds that of shortwave radiation ( 2 SW ) are included, yielding = 29 models. The gray dashed lines represent the best ODR fit between Δ a and a estimates across the models, with light gray shading indicating the 95 % prediction interval. The optimal temperature sensitivity corresponds to the point where the ODR fit line crosses Δ a = 0 eV (i.e., no adjustment to a is needed) and is indicated by a green dashed line. ODR, orthogonal distance regression.     Fig. 8 | Models for which NEE estimates explain atmospheric CO 2 variability less well than GPP estimates (red lines and symbols) show biased NEE seasonal cycles. Multimodel mean seasonal cycles of GPP (negative for uptake; dashed lines), ecosystem respiration ( E ; dotted lines), and NEE (negative for net uptake; solid lines) in (a) North America, (b) croplands, (c) evergreen needleleaf forests, and (d) deciduous broadleaf and mixed forests during 2007-2010 are shown for models for which the simulated GPP captures the variability in observed atmospheric CO 2 ( 2 GPP ) better than does shortwave radiation ( 2 SW = 0.23). Blue lines indicate seasonal cycles of carbon fluxes from terrestrial biosphere models (TBMs) from the MsTMIP v2 and TRENDY v6 ensembles for which NEE estimates outperform GPP estimates in explaining the observed CO 2 variability ( 2 NEE > 2 GPP ; = 16), whereas red lines represent seasonal cycles from TBMs for which GPP estimates outperform NEE estimates ( 2 NEE < 2 GPP ; = 10). Seasonal cycles of carbon fluxes from the FLUXCOM models ( = 3) are in orange. In addition, the geostatistical inverse model estimates of NEE (GIM NEE) are shown as a reference (gray). Months of peak fluxes are indicated by circles (NEE), upward triangles (GPP), and downward triangles ( E ). GPP, gross primary productivity; NEE, net ecosystem exchange.   Fig. 9 | Models for which NEE estimates explain atmospheric CO 2 variability less well than GPP estimates (red lines and symbols) show improved NEE seasonal cycles after rescaling ecosystem respiration estimates using the optimal temperature sensitivity for North America. Multi-model mean seasonal cycles of GPP (dashed lines), ecosystem respiration rescaled using the optimal temperature sensitivity for North America ( E ; dotted lines), and the resultant NEE (solid lines) in (a) North America, (b) croplands, (c) evergreen needleleaf forests, and (d) deciduous broadleaf and mixed forests during 2007-2010. Similar to Fig. 8, only models for which GPP estimates explain observed atmospheric CO 2 variability better than shortwave radiation ( 2 GPP > 2 SW = 0.23) are included. Blue lines indicate seasonal cycles of carbon fluxes from TBMs for which the explanatory power of original NEE estimates exceeds that of GPP estimates ( 2 NEE > 2 GPP ; = 16), whereas red lines represent seasonal cycles from TBMs for which NEE originally trails GPP in explanatory power ( 2 NEE < 2 GPP ; = 10). Seasonal cycles of carbon fluxes from the FLUXCOM models ( = 3) are in orange. In addition, GIM NEE is shown as a reference (gray). Months of peak fluxes are indicated by circles (NEE), upward triangles (GPP), and downward triangles ( E ). GPP, gross primary productivity; NEE, net ecosystem exchange.  Fig. 2), whereas horizontal error bars indicate ±1 standard deviation of mean annual air temperature (based on 328 grid cells for croplands, 544 for evergreen needleleaf forests, and 357 for deciduous broadleaf and mixed forests). Also shown is the relationship between climatological 10 for heterotrophic respiration derived from field observations and mean annual air temperature (gray line) from Koven et al. (2017) for comparison. Blue and red shaded areas indicate high-sensitivity, cold-climate emergent domain and low-sensitivity, warm-climate emergent domain (Koven et al., 2017), respectively.

Notes 1 Combination tests to assess impact of the uncertainty in gross primary productivity on the inferred temperature sensitivities of ecosystem respiration
To test the influence of uncertainty in model estimates of gross primary productivity (GPP) on the inferred temperature sensitivities of ecosystem respiration, we mix and match GPP and ecosystem respiration estimates drawn from different model simulations. For each "mix-and-match" pair we recalculate the adjustment to the activation energy of ecosystem respiration in the same way as described in the Methods in the main text (see "Optimizing the temperature sensitivity of respiration"). For quality control, the combination tests are limited to the 29 model simulations for which the explanatory power of GPP estimates ( 2 GPP ) exceeds that of shortwave radiation ( 2 SW ; Extended Data Fig. 6), consistent with Fig. 2 in the main text. This means that for each model simulation, the estimates of ecosystem respiration are paired with GPP estimates from each of the other 28 model simulations to obtain a set of optimal adjustments to the activation energy of ecosystem respiration under different GPP scenarios. Then, we obtain the average adjustment for each model simulation. Lastly, the North American and biome-specific optimal temperature sensitivities are derived from the same orthogonal distance regression method used to obtain the main results (Fig. 2 in the main text).
Results from the combination tests confirm that GPP uncertainty does not substantially impact the inferred optimal temperature sensitivities of ecosystem respiration. Differences between the optimal temperature sensitivities derived from the combination tests ( Supplementary Fig. 1) and those derived based on models' original GPP estimates (Fig. 2, main text) are well within the uncertainty range. The optimal temperature sensitivity for croplands derived from the combination tests (0.32 eV; Supplementary Fig. 1b) shows the largest difference from the original optimal estimate (0.38 eV). However, this is not surprising given that croplands are the largest contributor to model divergence in North American GPP estimates (Sun et al., 2021). Adjustments to estimates of activation energy for individual models in the combination tests are also highly consistent with adjustments based on each model's original GPP estimates ( Supplementary Fig. 2).

Notes 2 Potential impact of lateral emissions on the inferred temperature sensitivities
Terrestrial biosphere models (TBMs) vary in terms of their representation of carbon emissions from lateral transport processes (Huntzinger et al., 2013). Chief among these processes are consumption of crop products, decay of wood products, and outgassing from inland waters (Ciais et al., 2021). It remains difficult to incorporate these flux components into the evaluation of carbon budgets due to the lack of well-validated gridded observational data sets. Here, we construct a synthetic field of plausible emissions from lateral fluxes and evaluate the impact of these emissions on the inferred temperature sensitivities of respiration for North America and individual biomes. Because TBMs may resolve different components of lateral emissions and it is challenging to harmonize model differences, this test is intended to give an upper bound of the potential impact of lateral emissions.
We construct the synthetic lateral emission field based on the magnitudes of North American crop product (182 Tg C yr −1 ), wood product (76 Tg C yr −1 ), and riverine (219 Tg C yr −1 ) emissions compiled in Ciais et al. (2021) and the spatial distribution of crop, wood, and riverine lateral fluxes in Byrne et al. (2023). First, we assume that the lateral fluxes that represent additional efflux of carbon, i.e., carbon transported into a region from elsewhere (for example, crop and wood products shipped to population centers), are fully transformed into CO 2 emissions. Following this assumption, we use the spatial distribution of lateral fluxes from crop trade, wood trade, and riverine transport (Byrne et al., 2022;Byrne et al., 2023) to represent the spatial distribution of corresponding emissions. We then scale these emissions to the magnitudes given in Ciais et al. (2021). Due to the lack of monthly resolved emissions, we distribute lateral emissions equally across months. The resulting lateral emission field is then converted to atmospheric CO 2 enhancements using the same atmospheric transport footprints used for the main results (see "Optimizing the temperature sensitivity of respiration" in the Methods). We then subtract the influence of lateral emissions (mean: 0.27 ppm; median: 0.20 ppm) from atmospheric CO 2 observations. Finally, we perform the same optimization procedure to infer the optimal temperature sensitivities as we did for the main results (Fig. 2).
Inferred optimal temperature sensitivities for North America and individual biomes are indistinguishable from those derived without considering potential lateral emissions ( Supplementary  Fig. 3). The largest differences in the optimal temperature sensitivity estimate are found in forests (0.03 eV) but are still within the uncertainty range. Unsurprisingly, adjustments to estimates of activation energy for individual models in the lateral emission test also agree with adjustments without considering potential lateral emissions ( Supplementary Fig. 4).
These results suggest that different ways of accounting lateral emission (or lack thereof) in TBMs are unlikely to meaningfully affect inferred temperature sensitivities of ecosystem respiration at the biome scale.

Notes 3 Soil moisture influence on model-represented temperature sensitivities
To account for the potential influence of soil moisture on the temperature sensitivity of ecosystem respiration, we fit the following regression: where E (µmol m −2 s −1 ) is model-simulated ecosystem respiration, (m 3 m −3 ) is the mean volumetric soil water content in the top 2 m, obtained from North American Regional Reanalysis (Mesinger et al., 2006), 0 , 1 , and 2 are the coefficients, and is the residual. From the coefficient of 1/ we obtain the temperature sensitivity of ecosystem respiration, a,SM (eV), with soil moisture influences removed.
We find that accounting for soil moisture influences does not yield different estimates of modelrepresented temperature sensitivities of ecosystem respiration ( Supplementary Fig. 5). Although soil moisture modifies the temperature sensitivity of ecosystem respiration at the plot scale and on weekly to seasonal timescales (Reichstein et al., 2002), this effect is minor at the biome scale and over a longer time horizon. Thus, misidentification of the temperature sensitivity due to misrepresented soil moisture sensitivity seems unlikely.

Notes 4 Radiation influence on model-represented temperature sensitivities
Because shortwave radiation is highly correlated with temperature at a monthly resolution and on the biome scale, we take a two-step regression approach to account for the partial influence of shortwave radiation on model-represented temperature sensitivities of ecosystem respiration. To remove the co-variability between shortwave raidation and temperature, we first fit the following regression: where sw is the downward shortwave radiation flux at the surface (W m −2 ) obtained from North American Regional Reanalysis (Mesinger et al., 2006), 0 and 1 are regression coefficients, and sw is the residual.
We then fit the following regression to account for the partial influence of shortwave radiation in determining the temperature sensitivity: where 0 , 1 , and 2 are regression coefficients and is the residual. Similarly, from the coefficient of 1/ we obtain the temperature sensitivity of ecosystem respiration, a,SW (eV), with the partial influence of shortwave radiation removed.
We find that removing the partial influence of shortwave radiation does not meaningfully alter model-represented temperature sensitivities of ecosystem respiration on the biome scale (Supplementary Fig. 6). Although radiation is often examined as a covariate of ecosystem respiration (Jung et al., 2017), its causal link to ecosystem respiration is mediated by its effects on GPP (Monteith, 1972) and the correlation between GPP and ecosystem respiration (Baldocchi, 2008). Consequently, the absence of a substantial influence of GPP bias on the inferred temperature sensitivity of biome-scale ecosystem respiration ( Supplementary Fig. 1) indicates that it is unlikely that radiation would bias the inferred temperature sensitivities.

Notes 5 Influences of acclimation on model-represented temperature sensitivities
Plants and microbes have been found to physiologically adjust the temperature response of respiration according to ambient temperature, known as thermal acclimation (Atkin & Tjoelker, 2003;Bradford et al., 2008). To account for the impact of potential thermal acclimation on the temperature sensitivity of ecosystem respiration, we assume that the activation energy ( a , eV) is a linear function of the mean annual temperature: where Δ (eV K −1 ) is the linear temperature response of a , (K) is the mean annual temperature, a,0 (eV) is the temperature sensitivity at a reference mean annual temperature 0 = 273.15 K. This also means that a = Δ .
Substituting Eq. 4 in the Arrhenius equation for ecosystem respiration (Eq. 1 in the main text), we obtain where is the pre-exponential factor in the Arrhenius equation.
We then fit the regression From the coefficient of 1 we obtain the temperature sensitivity of ecosystem respiration at the reference mean annual temperature, a,0 . From the coefficient of 2 we obtain Δ , the temperature response of a .
Accounting for the thermal acclimation effect does not yield substantially different estimates of model-represented temperature sensitivities ( Supplementary Fig. 7). Therefore, thermal acclimation of respiration does not explain the spread in model-represented temperature sensitivity of ecosystem respiration, or the differences in such sensitivity between plot and biome scales and across biomes.

Notes 6 Analysis of the seasonal cycles of carbon fluxes
We find strong differences, both in terms of seasonality and magnitude, in estimates of net ecosystem exchange (NEE), gross primary productivity (GPP), and ecosystem respiration ( E ) between the group of TBMs for which NEE estimates explain atmospheric CO 2 variability better than those models' GPP estimates ( 2 NEE > 2 GPP ) and the group of TBMs for which moving from GPP to NEE degrades performance ( 2 NEE < 2 GPP ). The differences in the seasonal cycle of NEE between these two groups are especially notable in the croplands (Fig. 8b), showing up as a two-month early bias of the peak timing of NEE for TBMs for which GPP outperforms NEE ( 2 NEE < 2 GPP ; red lines and symbols) relative to that of the geostatistical inverse model estimates of NEE (GIM NEE; gray line and symbols; Shiga et al., 2018). This phase bias in cropland NEE is caused by both an offset of the timing of peak GPP and an overestimate of the seasonal amplitude of ecosystem respiration relative to TBMs for which NEE explains observed atmospheric variability better than GPP ( 2 NEE > 2 GPP ; blue lines and symbols) and FLUXCOM models (orange lines and symbols; Fig. 8b). By contrast, the seasonal cycles of NEE in deciduous broadleaf and mixed forests (Fig. 8d) do not show a phase bias but differ strongly in amplitude across model groups. Model estimates of NEE in evergreen needleleaf forests do not differ much among different model groups (Fig. 8c).
Analysis of the seasonal cycles of carbon fluxes further implicates the seasonal amplitude of ecosystem respiration estimates ( Fig. 8; dotted lines) as a source of the NEE phase bias that degrades performance. For TBMs for which GPP outperforms NEE in explaining atmospheric CO 2 variability ( 2 NEE < 2 GPP ), the one-month early bias in the peak timing of cropland GPP relative to the rest of TBMs and FLUXCOM models ( Fig. 8b; red dashed lines) is amplified to a two-month bias in the peaking timing of NEE (red solid lines) by the high bias in the seasonal amplitude of ecosystem respiration (red dotted lines). This is also true for the North American domain (Fig. 8a), where the high bias in the seasonal amplitude of ecosystem respiration causes a phase bias in NEE for TBMs for which GPP outperforms NEE ( 2 NEE < 2 GPP ), despite the fact that the seasonal cycles of GPP are nearly identical between the two groups of models. In forested biomes (Fig. 8c, d), there is also a consistent high bias in the seasonal amplitude of ecosystem respiration for TBMs for which GPP outperforms NEE ( 2 NEE < 2 GPP ), even though this bias does not cause a phase bias in NEE. Moreover, for these TBMs ( 2 NEE < 2 GPP ), the high bias in peak respiration is always accompanied by a smaller but noticeable low bias in winter respiration relative to TBMs for which NEE outperforms GPP ( 2 NEE > 2 GPP ) and relative to FLUXCOM models (Fig. 8). This suggests that the amplitude bias in ecosystem respiration stems from the temperature sensitivities of ecosystem respiration as represented within each model.
Rescaling ecosystem respiration by the ensemble-optimal temperature sensitivity over North America (0.43 eV) leads to reduced bias in both the seasonal amplitude of ecosystem respiration and the timing of peak NEE for models for which NEE performance trailed that of GPP ( 2 NEE < 2 GPP ; red lines and symbols) in the North American domain (Fig. 9a) and croplands (Fig. 9b). For these models, the NEE phase bias relative to GIM NEE is eliminated over North America and reduced to one month in croplands after rescaling. In forested biomes (Fig. 9c, d), however, rescaling causes inconsistent shifts of the NEE peak from June to July. But this is mainly due to the small difference in NEE between June and July in forested biomes. Overall, the comparison of the seasonal cycles of NEE and ecosystem respiration before and after the rescaling lends strong support to the finding that bias in the temperature sensitivity of ecosystem respiration is a leading cause of NEE underperformance.