Numerical simulations

Figure 1 presents implications of new R d components of Eq. (1). Figure 1a shows for broadleaf trees significant increases across all temperatures in respiration compared to standard JULES, when using the new R d,25 values (T G = 25 °C) and either Q 10 = 2 or the b,c T l response. Figure 1b shows the four PFT responses to T l , with revised R d,25 values, T G again 25 °C, and b,c formulation. Figure 1c illustrates strong R d differences of Eq. (1) between acclimation temperatures T G = 15, 25 and 35 °C (for broadleaf trees). In Fig. 1d, the orange curve is the same R d −T l response (T G = 25 °C) as in Fig. 1c. However, the red curve sets acclimation temperature equal to instantaneous temperature i.e. T G = T l . This sensitivity test recognises that although acclimation growth temperature, T G , is determined over longer 10 day periods, higher T G values will be geographically where T l is higher and vice versa. This dampens R d variation in T l . Curve dashed for extremely rare temperatures T G > 35 °C.

Fig. 1 Upper canopy leaf-level dark respiration. a Standard Joint UK Land Environmental Simulator (JULES) model (with TRY-based σ l and n l0 values and Q 10 response modulated at high and low temperatures (Methods)). Also revised R d,25 (T G = 25 °C) with both Q 10 = 2.0 and b,c responses to T l . b New R d,25 and b,c response to T l , for other PFTs. C 3 grasses near identical curve to Shrubs. c New R d,25 and b,c formulation for broadleaf trees, but alternative acclimation temperatures. d New R d,25 and T l as b,c again broadleaf trees, for both T G = 25 °C and T G = T l . Orange curve common all panels. Light inhibition not included in responses Full size image

JULES scales R d to canopy-level respiration, R d,c (μmol CO 2 m−2 s−1). It can calculate CO 2 exchange at each canopy level19, including dependence on vertical decline of leaf nitrogen19 and differentiation of direct and diffuse radiation20. However, data are unavailable for how well Eq. (1) performs at lower canopy levels, even if nitrogen concentration and temperatures are known. Given this, we use a simple big-leaf exponential decline in leaf respiration throughout the canopy, decay co-efficient k = 0.5 and dependent on leaf area index (LAI). Implicit is that canopy nitrogen and light levels decay identically, affecting respiration and photosynthesis. A 30% light inhibition of non-photorespiratory mitochondrial CO 2 release21 is included for light above ~ 2 W m−2. R d,c is also reduced by any soil moisture stress4. Other respiratory components4 include maintenance respiration of stems and roots. These are modelled as larger for canopies with high LAI, higher nitrogen concentrations and also described as scaled in R d,c . Combining stem and root respiration with R d,c gives whole-plant maintenance respiration, R pm (Methods). JULES calculates growth respiration R pg as linearly increasing (via co-efficient r g ) with plant GPP minus R pm . Combining R pm and R pg gives whole-plant respiration, R p (μmol CO 2 m−2 s−1). Hence changes to R d description influence all respiration components that add together to give R p .

Spatial gridded calculations enable geographical implications of revised R d description to be determined. JULES is first forced by meteorological conditions near to pre-industrial state, using UK Climate Research Unit data and atmospheric CO 2 concentration of 280 ppm. To understand the implications of each part of our new R d description and as given by Eqn. (1), we add sequentially each component to understand its relative importance. Four sets of simulations are as follows. The first is called Standard—this is standard default JULES (although with TRY-based n l0 and σ l values, as in all simulations; Q 10 response but high/low temperature suppression—see Methods). The second is called New_R d,25 —this is new baseline alone, hence new R d,25 values5, plus Q 10 = 2.0 (no suppression) T l response and fixed T G = 25 °C. The third is called New_R d,25 _b,c —this is new baseline and instantaneous response i.e. new R d,25 values and T l response with b,c formulation, but still T G = 25 °C. The fourth is called New_R d,25 _b,c_acclim—this is all factors including acclimation, i.e., new R d,25 values, b,c formulation and acclimation via variation in T G . The fourth simulations are therefore for the full Eq. (1)5, 6. Figure 2a–c shows how each of these new components of our respiration function uniquely influences whole-plant respiration globally for pre-industrial climate forcings. Figure 2a shows annual gridbox mean R p (weighting PFTs by fractional cover) for New_R d,25 minus Standard simulations. This shows that altered baseline through the new R d,25 (and removed suppression) causes large R p increases, including at mid-Northern latitudes and especially for the tropics. Figure 2b shows New_R d,25 _b,c minus New_R d,25 , illustrating the implications of the new instantaneous temperature description. The b,c formulation suppresses respiration in mid-latitudes but enhances for the tropics, although changes are smaller than Fig. 2a. Figure 2c presents New_R d,25 _b,c_acclim minus New_R d,25 _b,c, showing acclimation introduction generally increases predicted pre-industrial R p , except in the tropics where acclimation to higher temperatures lowers respiration.

Fig. 2 Gridbox-mean maps of total plant respiration for new processes and imposed climate change. Changes to R p as: a introduction of new R d,25 values: New_R d,25 minus Standard; b effect of new instantaneous T l response: New_R d,25 _b,c minus New_R d,25 ; c effect of acclimation: New_R d,25 _b,c_acclim minus New_R d,25 _b,c; d effect of climate change to present, as year 2015 minus year 1860 and new processes New_R d,25 _b,c_acclim; e similar to d, but year 2050 minus year 2015. Panel f actual 2015 values of New_R d,25 _b,c_acclim. Scales different between panels to highlight effects. Units are SI. Panels d–f are means across the 34 GCMs emulated Full size image

To estimate anthropogenically-induced climate change, changes in near-surface meteorological conditions use the Integrated Model Of Global Effects of climatic aNomalies (IMOGEN) pattern-scaling system22, 23 (Methods) responding to altered atmospheric greenhouse gas (GHG) concentrations. Patterns are calibrated against 34 GCMs of the Coupled Model Intercomparison Project Phase 5 (CMIP5)24, with identical method to that originally undertaken for the HadCM3 GCM22. We use known historical, then future GHG concentrations of the RCP8.5 business-as-usual representative concentration pathway25. Figure 2d shows, for the most complete updated model, New_R d,25 _b,c_acclim, that historical climate change increases R p in most locations and especially tropics, despite acclimation dampening stimulatory warming effects. Figure 2e presents calculations between years 2015 and 2050, showing R d with similar changes to recent past in Fig. 2d. Figure 2f presents year 2015 absolute R p values for New_R d,25 _b,c_acclim case.

Figure 3 presents model output for single illustrative locations and in year 2015. For our four simulations, presented are respiration components (R d , R d,c and R p ), plus GPP and NPP. We chose seven sites across South America, a temperate grassland (London) and boreal region shrubs (Siberia). We select multiple South America sites (Methods), as these are some of the few where measurements are available of all respiration components. In general, new R d,25 values, whether also with or without adjustment by the b,c formulation and acclimation, give marked increases in predicted respiration. Transition to whole canopy (R d,c ) and whole plant (R p ) respiration illustrates how our leaf level changes propagate to these aggregated fluxes. Uncertainty bounds5 on r 0 , r 1 and r 2 are propagated through the JULES model (Methods) to give uncertainty on R d,c and R p as shown in Fig. 3, while measurement uncertainty is from the literature describing each site. For South American sites, and with our choice of big-leaf approximation, our changes reproduce whole-canopy respiration R d,c better (i.e., model and data uncertainty bounds overlap, and better than the default Standard JULES configuration), and in some instances also R p . More specifically, we define the JULES model as having improved performance when the Standard simulation estimate of R p lies outside the data-based bounds on whole-plant respiration, but simulations New_R d,25 _b,c_acclim then fall within those bounds. This happens for the sites at Manaus, Tambopata, Iquitos (dataset a), and Guarayos (dataset a). However, when subtracting R p from GPP estimates, NPP values are generally too small. We note that observations of nitrogen at different canopy positions from tropical tree species suggest an effective decay co-efficient k with value nearer to 0.2 than 0.526. Using this to scale, and with Eqn (1) still used for upper canopy levels, gives exceptionally large R p values and unsustainable negative NPP.

Fig. 3 Respiration and primary productivities (gross as GPP and net as NPP) at selected locations during modelled year 2015. Seven locations (details in Methods) a are South American b–h, along with i for gridbox containing London, UK, and j is in Siberia, Russia (Lat 70 N, Lon 82.5 E). Shown for dominant plant functional type (PFT) at each site, left to right, for each histogram cluster: upper canopy leaf-level respiration (with light inhibition) R d , whole canopy-level respiration R d,c , total plant respiration R p , GPP and NPP. Each histogram cluster are four estimates: Standard, New_R d,25 , New_R d,25 _b,c and New_R d,25 _b,c_acclim. South America sites, 5th (or 6th) column are measurements. Dominant PFTs: b–h BTs, i grasses, j shrubs. Uncertainty bounds of ± one s.d. are presented which for model estimates are from data-based upper canopy leaf-level uncertainty estimates, subsequently propagated through the model. For measurements, these bounds are taken from the literature Full size image

Figure 4 shows global time-evolving changes, since pre-industrial times, in total whole-plant respiration, ΔR p (GtC yr−1) and for our four RCP8.5 scenario simulations. Annotated are pre-industrial and 2015 absolute R p estimates. Replacement of standard JULES with GlobResp-based5 R d,25 values (still Q 10 = 2, although with no high or low temperature suppression) approximately doubles both pre-industrial respiration estimates (as marked) and the projected changes in ΔR p under climate change. Replacing Q 10 with b,c formulation6 causes slight global changes. Thermal acclimation increases R p slightly for pre-industrial but decreases evolving ΔR p , i.e., comparing simulations New_R d,25 _b,c_acclim and New_R d,25 _b,c. The stimulatory effect of acclimation arises from the higher predicted rates in globally widespread biomes where T g < 25 °C, but then dampens responses of such sites to future warming. Our new global R p values (80.4 GtC yr−1 in 2015 for New_R d,25 _b,c_acclim simulations) are higher than other estimates for contemporary periods. One27 global GPP estimate is 119.6 GtC yr−1, and balancing soil plus plant respirations will be similar magnitude i.e. together they are also of order 120 GtC yr−1. With soil respiration equivalent in size to R p , this suggests plant respiration of order 60 GtC yr−1. Our analysis implies global R p values could be ~30% higher than that value. However, this estimate is not just a consequence of the entrainment of GlobResp data in to the JULES model, but also the scalings within it and as illustrated at selected geographical points in Fig. 3.

Fig. 4 Time series of change in areally-averaged global respiration. Presented are time-evolving model estimates of change in total whole-plant respiration, ΔR p . The colours of turquoise, blue, yellow and magenta are Standard, New_R d,25 , New_R d,25 _b,c and New_R d,25 _b,c_acclim respectively. Where yellow and blue projections overlap, the colour is brown. The spread corresponds to the different projections of climate drivers, based on the 34 Global Circulation Models (GCMs) emulated in the Integrated Model Of Global Effects of climatic aNomalies (IMOGEN) modelling system and for RCP8.5 scenario. The continuous lines are the mean, and the spread as ± two s.d. which broadly covers inter-GCM spread. Pre-industrial (marked Pre-I) and year 2015 model mean absolute estimates of R p are as annotations Full size image

The GlobResp database5 is sufficiently comprehensive as to be globally representative for our simulations. Our analysis has implications for other land ecosystem modelling groups. From a survey of ten leading land surface models, six of these simulate leaf respiration with a dependency on nitrogen content (models listed in Methods). In addition as is common to most land surface carbon cycle models (e.g., also for the Lund Potsdam Jena LPJ model (Table 1 of ref. 28)), the JULES system calculates maintenance respiration for other components of roots and stems that are based on their carbon and estimated nitrogen content. This approach follows pioneering research29 which moved respiration representation on from simply a fraction of GPP, as had been assumed beforehand. We expect the impact of changing the temperature function itself to at least be common among the current generation of models. However this does open research questions as to how Eq. (1) might change at lower positions in canopies, and whether root, stem and growth respiration models require refinement. This is especially because our GlobResp-based changes propagate directly through modelled canopies by JULES (Eqs. 42–25 in ref. 4). Hence higher upper-canopy R d values then generate larger rates of whole-plant respiration, R p , than other estimates.

Benchmarking tests of modelled respiration fluxes will be important. For instance, the International LAnd Model Benchmarking project (ILAMB)30 is a comprehensive system collating datasets relevant to land surface functioning and of importance to land surface respiration is the Global Bio-Atmosphere Flux (GBAF)31 dataset based on extrapolation of eddy-covariance FLUXNET sites. Also available are estimates of global soil respiration32, which in conjunction with GBAF measurements can return total plant respiration, at least for comparison at night-time periods. Presently, however, without comprehensive measurements of other canopy components, it is difficult to attribute any discrepancies to GlobResp versus lower-canopy, stem, root or growth contributions. Should higher R p values imply especially low values of NPP, then GPP parameterisation may need reassessment; other analyses suggest current estimates of GPP may be too low33. Despite this, in Fig. 5 we perform large-scale comparisons against two Earth Observation-based datasets. These are estimates of NPP from the MODerate-resolution Imaging Spectroradiometer (MODIS) satellite, using the MOD17 algorithm34, 35, and of GPP from the Model Tree Ensemble (MTE) method27. For both datasets, we evaluate mean NPP and GPP values depending on location, and mapping these on to local dominant biomes based on the World Wildlife Fund (WWF) ecoregion classifications36 (Methods). These data-based estimates locally represent mean NPP and GPP, and so for parity we compare against modelled gridbox mean JULES calculations of the equivalent quantities. That is, we use areal weighting of the five PFT types in JULES for each position. To keep similarities with the WWF categories, we plot in Fig. 5 total annual NPP and GPP for both data and JULES, integrated over areas for the named biomes as marked. Presented are Standard and New_R d,25 _b,c_acclim simulations. Calculations with New_R d,25 and New_R d,25 _b,c model format are very similar to New_R d,25 _b,c_acclim and so not shown. As expected, in all cases, introduction of GlobResp-based respiration estimates results in much lower modelled NPP values. Furthermore for New_R d,25 _b,c_acclim simulations and all eight biomes, these are significantly less than MODIS-based measurements. The two sets of simulations have similar GPP estimates, illustrating weak indirect couplings in the JULES model between respiration changes and influence (e.g., via hydrological cycle) on gross primary productivity. It is noted in Fig. 5b that JULES model estimates of GPP are similar to the MTE-based data for tropical forests and tropical savannahs. Uncertainty bounds on data adopt the global literature values of ± 15% for NPP37 and ± 7% for GPP38. These are the small horizontal black bars, shown only on New_R d,25 _b,c_acclim red points.

Fig. 5 Data- and model-based global estimates of net primary productivity and gross primary productivity for different biomes. a Global measurements of total annual mean net primary productivity (NPP), average for years 2000–2011, and using Earth-observed MODerate-resolution Imaging Spectroradiometer (MODIS) measurements. Values are spatially aggregated for different World Wildlife Fund (WWF) biome classifications. The dominant biome type at each location is linked to NPP with the MOD17 algorithm applied to MODIS values (horizontal axis). Similarly gridbox-mean JULES estimates of NPP are multiplied by gridbox area, and combined for each WWF biome (vertical axis). This is dependent upon which WWF biome is dominant for the grid location. Note logarithmic axes. JULES NPP estimates are slightly negative for Mediterranean grasslands and so off axes. b Similar calculation for gross primary productivity (GPP), with measurements from the Model Tree Ensemble (MTE) algorithm. Both panels, model values presented in blue for standard JULES version (i.e., Standard simulation), and in red for new R d,25 values with b,c temperature response and acclimation (i.e., New_R d,25 _b,c_acclim simulation). For GPP, differences are small between two model forms, with red symbols overlapping blue symbols. Uncertainty bounds on NPP and GPP data are the small black horizontal bars ( ± one s.d.), shown for red symbols only. All calculations include only land points with less than 50% agriculture Full size image

In Fig. 6, we add geographical information to our global data estimates of NPP and GPP, and for corresponding JULES simulations with all effects, i.e., New_R d,25 _b,c_acclim (expanding on the red symbols of Fig. 5). Figure 6a is JULES NPP estimates divided by MOD17-based NPP estimates (and multiplied by 100 to give percentage). In general modelled NPP with new plant respiration description, is smaller than MOD17 NPP across the geographical points. For some points it can give unsustainable negative modelled NPP values. For GPP, the situation is slightly less clear. Figure 6b is JULES GPP estimates divided by MTE-based GPP values, again as percentage. For many points, the JULES model is also underestimating GPP, and this includes much of the Amazon region. However, for the tropics, a few modelled GPP values are actually higher than data. This offers an explanation as to why GPP appears underestimated in some tropical points of Fig. 3, yet for the average across Tropical Forest (TF), JULES performs well (Fig. 5b). Figure 6b also shows that modelled GPP is usually too low outside of the tropics. This is why, when combined with the enhanced respiration of New_R d,25 _b,c_acclim formulation, this can lead to very low or even unsustainable negative NPP. Figure 6c shows the dominant WWF-defined biomes for each location.