Updated respiration routines alter spatio-temporal patterns of carbon cycling in a global land surface model

We updated the routines used to estimate leaf maintenance respiration (MR) in the Energy Land Model (ELM) using a comprehensive global respiration data base. The updated algorithm includes a temperature acclimating base rate, an updated instantaneous temperature response, and new plant functional type specific parameters. The updated MR algorithm resulted in a very large increase in global MR of 16.1 Pg (38%), but the signal was not geographically uniform. The increase was concentrated in the tropics and humid warm-temperate forests. The increase in MR led to large but proportionally smaller decreases in global net primary production (19%) and in average global leaf area index (15%). The effect on global gross primary production (GPP) was a more modest 5.7 Pg (4%). A detailed site level analysis also demonstrated a wide range of effects the updated algorithm can have on the seasonal cycle of GPP. Output from the updated and old models did not differ markedly in how closely they matched a suite of benchmarks. Given the substantial impact on the land surface carbon cycle, a neutral influence on model benchmarks, and better alignment with empirical evidence, an MR algorithm similar to the one presented here should be adopted into ELM.


Introduction
A substantial portion of the Earth's mass and energy balance is controlled by land surface vegetation through links in the energy, water, carbon and nutrient cycles. For example, total land-based autotrophic dark respiration is estimated to release around 60 Pg of carbon each year [1][2][3]. A major component of autotrophic dark respiration is respiration taking place in mature, fully expanded leaves (R d ), with most of the energy produced by R d being used to support cellular maintenance processes in leaves. However, a mechanistic model of dark respiration has proved elusive. Thus, algorithms in land surface models used to simulate the global carbon budget must rely on empirical models of respiration, which are in turn dependent on the data used to derive them. Here, we show that combining the two most comprehensive analyses of global respiration to date in an updated leaf maintenance respiration (MR) algorithm has a substantial influence on the global estimate of maintenance respiration, which propagates through the entire global carbon cycle.
The importance of respiration and lack of consensus on the form that algorithms can take has been well documented-for example [4], reported 15 identifiable approaches from 21 different LSMs. A common approach is to use respiration-nitrogen scaling to estimate a base rate of leaf R d at a standardized temperature. However, a fixed respiration-nitrogen relationship misses large variations in base respiration across broad categories of plants for example, R d has been reported higher in forbs and grasses than woody angiosperms for a given foliar nitrogen concentration [5,6]. Thereafter, leaf R d is often modeled to vary with temperature assuming a fixed Q 10 [7]. The fixed Q 10 response ignores the importance of temperature acclimation on both long (weeks to months) and short (minutes to hours) timescales [8,9], The importance of acclimation has been highlighted, but not yet incorporated into the primary code of LSMs [7,[10][11][12][13][14]. The adoption of acclimating algorithms has likely been resisted, in part, due to the interlinked nature of LSM code and the need to modify other components, e.g. [15]. Here, we show that respiration with both long and short term temperature acclimation can be incorporated into a land surface model (the Energy Land Model, ELM) with negligible adjustments to other model components.
The current base respiration algorithm in ELM is based on nitrogen concentration data from a small pool of 11 species [16]. In turn the temperature response is identical to that of photosynthesis [17]. There is now a large body of empirical evidence for variation in R d that goes well beyond these approaches. For example [18], described a 16-fold range in leaf R d at a common temperature from 20 sites around the globe. Additionally, there appears to be nearly as much variation in leaf R d among species that co-occur within sites as across biomes, implying both genetic and environmental controls on leaf R d [19][20][21][22][23][24][25]. Acclimation to temperature has also been documented in a number of studies [5,23,26]. Several recent physiological advances [5,8,9,27] have formally quantified these sources of variation in a manner that is amenable to incorporation in an LSM. These advances are timely, as the modeling community has also recognized that plant R d variability is not sufficiently captured in global land models [28][29][30][31].
Here, we focus on incorporating into ELM the respiration variation reported in two global studies. In [5] a dataset of global leaf R d was compiled. The global respiration (GlobResp) database consists of information from 899 species sampled from 100 sites around the world. This is, to the best of our knowledge, the largest data set of its kind. Thereafter [8], used data from 231 species distributed across seven biomes to demonstrate that leaf R d systematically exhibits a rapid temperature response that differs from the common fixed Q 10 . Further, the response is remarkably consistent across biomes and PFTs. Together, these studies represent the best available evidence of how leaf R d varies across organisms and responds to temperature.
A previous study used the JULES model to incorporate these algorithms [11], but was forced to modify the canopy nitrogen profile to stabilize net primary production (NPP). Such modifications may be necessary, but confound the interpretation of changes to model output from the new algorithms. Here, we make no modifications to transient model parameters other than those described in the updated leaf MR algorithms. This is particularly critical to this analysis, which focuses on the implications of the updated algorithms to multiple components of the land surface carbon cycle.
Given that the new findings highlighted above are based on almost two orders of magnitude greater data than the data base used in the current respiration algorithm in ELM, and capture biological processes that are seen as critical by the modeling, experimental, and field research biological communities, the time is right for these advances to be incorporated into a land surface model. After updating the ELM respiration algorithm we conduct a suite of analyses to show how the influence of MR runs through the global carbon cycle. First, we compare the global spatial patterns of the modified model to the default model; next we evaluate seasonal variation of the modified model relative to the default model and empirical observations at ten flux tower sites; and, last we compare both the modified and default models against global benchmarking data.

Model description
The U.S. Department of Energy's Energy Exascale Earth System Model (E3SM) is a fully coupled ESM with atmosphere, land and energy, ocean, sea ice, and land ice components. The land component of E3SM (ELM v1) is based on the Community Land Model 4.5 (CLM4.5) [32], and includes a canopy integrated leaf photosynthesis module [15,33], coupled carbon and nitrogen dynamics for vegetation and soil [34], vertical resolution of soil biogeochemistry [35] and a permafrost hydrology component [36]. In addition, it uses an improved photosynthesis module with temperature revised Rubisco kinetics from [37,38] the high temperature photosynthetic stress modifications from [17] and a prognostic phosphorus cycle [39,40]. The ELM model was run using fixed surface meteorology forcing data from GSWP3 [41] with 20 years of meteorology data (1901-1920), repeated over a 600 year spin-up followed by a 250 year accelerated spin-up [35,42] to stabilize carbon and nutrient pools. A transient simulation with historically varying atmospheric CO 2 concentrations, landuse change, and nitrogen deposition inputs starting at 1850, with the GSWP3 meteorology from 1901 to 2010, and detailed output from 1961 to 2010. For site-specific simulations, we used the global spin-up for the respective grid-cell, and climate data obtained from each flux site for the transient simulations. The model was seeded with approximately double the default initialization carbon to prevent depletion in tropical regions during model spin-up.

Respiration algorithms
The base rate in ELM's current respiration algorithm is based on empirical leaf R d -nitrogen relationships reported in [16] then modified by [43]. The current temperature response was first formulated in [17] and updated in [33]. These are combined to create the following mathematical form: where , nl pft is a plant functional type (PFT) grid cell specific massbased leaf nitrogen concentration (gN leaf m −2 ) which linearly modulates the base rate, C normalizes basal respiration to 25 • C, H a is activation energy (J mol −1 ), R is the gas constant, T 0 is 298.2 K, T is leaf temperature in kelvin, S v is the entropy value (J mol −1 K −1 ), and H d is deactivation energy (J mol −1 ). For clarity we update the algorithm in two steps, in accordance with the two studies that update the base rate of respiration and the instantaneous temperature response. To evaluate the effects of PFT-specific base rates and leaf-level R d -nitrogen relationships as well as incorporating medium-term temperature acclimation we used the GlobResp data, described in [5]: where r base,acclim is the base rate of leaf R d at a standard temperature of 25 • C, r 0,pft is the PFT specific coefficients related to the base rate, while r 1 and r 2 are the uniform nitrogen response and medium term temperature acclimation respectively. The nl pft term is mass based foliar nitrogen and T g ( • C) is the preceding 10 days average growth temperature, both of which vary by grid cell. Next, we updated the instantaneous leaf R dtemperature with an exponential function derived from [8]: Here, R d is the MR analogous to that from equation (1), r base,acclim is the acclimating base rate from equation (2) and T is the leaf temperature, a and b are empirical constants equal to 0.1012 and 0.0005. From here out, equation (3) is referred to as the variablebase-exponential (VBE) algorithm [44]. See table 1 for PFT specific parameter values and units. Leaf MR is upscaled to the canopy based on canopy leaf area index (LAI) and average nitrogen content. Recent work [45] has updated these routines, but the essential feature is an exponential relationship with LAI: In equation (4) n a is a dimensionless coefficient to evaluate canopy average nitrogen content, k n is an extinction coefficient and LAI is canopy LAI (m 2 m −2 ). The updated MR (equation (3)) is combined with the canopy scaler (equation (4)) and the total canopy LAI to provide an estimate of the canopy level maintenance respiration, R d,can : Further modifications to maintenance respiration, such as those due to other environmental factors like soil available water, are analogous to changes to the maximum carboxylation velocity (V cmax ). See [32] for complete details. Note that leaf MR is only one component of land surface respiration. Other living tissues also require maintenance respiration, which in ELM include: live stem, coarse roots, and fine roots. We do not update the respiration algorithms for these components and they operate under the following algorithm: 20 10 Q10 .
Here, R m,x (gC m −2 s −1 ) is the MR for a component x, which is either living stem, coarse roots or fine roots. They share a common base respiration r base (gC gN −1 m −2 s −1 , analogous to equation (1)) and each compartment has a unique nitrogen pool N x (gN m −2 ). They have a Q 10 temperature response determined by a local temperature for each of the pools, T x ( • C). The four land surface components with MR share a common growth respiration (GR) response, which is a fixed 30% of the total carbon allocated to each pool. The changes made to leaf MR have effects that ripple through these components despite no direct modification.

Benchmarking
We compiled eddy-covariance data from ten geographically and ecologically diverse FLUXNET sites. These sites were selected to represent a range of

Respiration algorithm formulations
The influence of the updated VBE algorithm on leaf R d is summarized in a set of idealized scenarios presented in figure 1. We fixed the foliar nitrogen content at 1%, 2%, and 3% (dry mass), corresponding to figures 1(a)-(c), and examined the response across a range of instantaneous temperature values. The updated algorithm and the default algorithm vary considerably across all nitrogen and temperature values, but the differences are magnified at lower nitrogen content and hotter temperatures. When the updated algorithm is aggregated across PFTs and held at a constant growth temperature of 25 • C, the simplest comparison to the current default, the updated estimates of respiration are doubled near 0 • C. As the temperature increases to 30 • C this difference magnifies to three times higher at 3% leaf N content (figure 1(c)) and six times higher at low leaf N content ( figure 1(a)). The PFT specific differences between the algorithms were greatest for grasses and lowest for needleleaf trees-as one would expect according to their respective base respiration rates (table 1). , approximately half the magnitude and proportion of change in MR. The NPP decreases were largest in the tropics but also occurred in the relatively colder and drier regions of both northern and southern hemispheres (figure S1). As expected, given the decrease in carbon available for NPP the global average LAI decreased by a comparable amount (15%) over the same time period (figure S2). The largest decreases were in the tropics, with smaller decreases in high latitude regions. However, despite decreases in NPP in middle latitudes there were some modest LAI increases there (figure S2).

Global analysis
The updated (VBE) model lowered global GPP by 5.7 Pg C yr −1 (figure 3), which was an order of magnitude lower effect (4% decrease) relative to MR (38% increase). Alterations in LAI patterns appeared to be the principal driver of GPP patterns, as they are closely aligned (figure S2). There are some modest increases in GPP across mid latitudes, but these are much smaller than the declines in GPP across the tropics and high latitude regions ( figure 3).
The overall modifications on the ELM estimated carbon cycle resulting from the VBE algorithm are summarized in figure 4. The increase in MR drives a smaller pool of carbon available for allocation, this in turn reduces the GR costs, which is why the decline in NPP is smaller than the decline in MR. However, the decrease in NPP is nearly proportional to the decline in LAI. Despite the substantial decline in global LAI, the concentration of this decrease in the tropics with generally high LAI values means that the majority of the decrease is concentrated in shaded LAI, with less influence on GPP than sunlit LAI. Note that the small global totals for sunlit LAI are a consequence of averaging over the entire globe over long time periods thus incorporating night time and regions with very low total LAI. In aggregate these changes cascade to the more modest decline in global GPP relative to the substantial increase in MR brought about by the VBE Figure 1. Foliar nitrogen effects on base respiration-instantaneous temperature relationships for default ELM respiration (e3sm-all veg), aggregated vegetation PFT using variable-base-exponential (VBE) parameters (vbe-all veg), and PFT specific VBE parameters (broad leaf trees, grasses, needle leaf trees, and shrubs). Foliar nitrogen (as a percent of leaf dry mass) levels for panels are 1% (a), 2% (b), and 3% (c) respectively. Aggregate PFT relationships are described by the dash-dotted lines (ELM-orange, VBE-brown) while PFT specific relationships for the VBE algorithm are shown in solid lines. The changes to above ground biomass are also reflected in a decline in the ratio of NPP to GPP [47], figures S3 and S4. This suggests that under the VBE algorithm there is less carbon that remains within global vegetation relative to the total amount fixed through photosynthesis. These declines are most severe in low productivity desert regions where the VBE algorithm reduces already marginal NPP values to nearly zero and by marked declines across tropical regions and eastern Siberia in the boreal zone ( figure  S4). The NPP/GPP ratio of the model is at the lower end of other global estimates [48,49], but within the boundaries of empirical studies [50,51]. The tropics and eastern Siberia also show notable declines in soil carbon stores (figure S5), which are most pronounced in the Arctic.

Flux sites
A detailed analysis at ten flux tower sites in varying climates found that the impact of the VBE algorithm was modest when averaged over the whole time period and growing season (figure S6). As expected from the global analysis the VBE algorithm had neutral or slightly lower annual GPP at most sites but at two temperate mixed forest sites (US-UMB and FR-Pue) annual GPP is greater ( figure S6). The modest changes in annual average GPP are contrasted by considerable variation over the course of the growing season, which is strongly dependent on the environment ( figure 5). Sites fell into roughly four different categories and an emblematic member of each  category is shown in figure 5: (a) GPP from the VBE algorithm was lower than from the default algorithm in all growing months; as found in tropical, temperate and boreal sites: BR-Sa1 ( figure 5(a)), GF-Guy, US-PFa, US-Ha1, and FI-Hyy; (b) at two warm temperate mixed forest sites, Fr-Pue ( figure 5(b)) and US-UMB, GPP using the VBE algorithm was higher than from the default at the peak of the growing season (June-September); (c) no change was found at the South African site; (d) the two Russian arctic sites, RU-Cok (figure 5(d)) and RU-Fyo showed subtle differences over the course of the short growing season with the VBE algorithm increasing GPP in the early growing season but decreasing it towards the end. See figure S7 for the seasonal cycle of sites not shown in figure 5. The differences between the algorithms are subtle enough that neither is a better match to the flux tower GPP, and generally the model shows less inter-annual variability than the flux tower estimates.

Benchmarking
Relative to default ELM respiration equations, the VBE algorithm had little to modest impact on most global ILAMB mean state benchmarking scores (table S2). Although differences were modest, the VBE algorithm more closely matched benchmarking scores than did default formulations for Biomass, CO 2 , and LAI while the old algorithm more closely matched Global Net Ecosystem Carbon Balance and Soil Carbon. However, as is clear from the global maps (figures 2, 3 and S1-S5) even in cases where VBE and default algorithms are comparable for a global aggregated benchmark, there are regional differences. For example, when comparing algorithms against the Global Biosphere-Atmosphere Flux GPP benchmark the VBE respiration algorithm reduced GPP over large portions of Brazil relative to the default, bringing portions of the upper Amazon basin much closer to the benchmark data, but performing worse across sections of the lower Amazon.

Summary
The changes implemented to the leaf MR algorithms in ELM bring the model into alignment with updated global empirical data. The data bases that informed the updated algorithms comprise a nearly two orders of magnitude increase in both the number of species sampled and the geographic coverage. The processes that we incorporated include a PFT dependent base respiration rate, medium time-scale temperature acclimation, an updated nitrogen-respiration relationship and a new instantaneous temperature response. When similar adjustments were made to the JULES model [11], a comparable change in global MR was found, but geographical patterns differed and changes to the canopy nitrogen profile were required to enable the model to run and to stabilize NPP. Related work on the Community Land Model focused on changes to respiration in conjunction with updated photosynthesis [13] and did not have access to the most recent MR databases that we used here. The improvements presented here bring the respiration algorithm in a global LSM up to the state of the art in empirical estimates.

Conclusion
The improvement in model algorithm fidelity to state-of-the-art empirical data had a dramatic effect on global autotrophic MR but surprisingly little influence on how well the model compared to a range of benchmarks, which themselves contain considerable uncertainty. Moreover, though our updated model and the prior default model did not vary markedly from one another in terms of annual averages across the globe or at any of the sites, we found considerable differences spatially around the globe (figures 2, 3, and S1-S5) and temporally across the growing season (figure 5). One might view the modest change in the annual average as encouraging, in that a substantial change to just one aspect of the carbon cycle of the model did not degrade model performance (as noted in [52]), but it is surprising that there was little improvement in matching the model output to benchmarks, given large advances in respiration parameterization and routines. If plant R d components of land surface models historically underestimated maintenance respiration, then other components of such models may have been tuned such that overall carbon cycle outputs better matched benchmark data. However, we would have expected greater degradation of model output if this were the case, as shown previously for boreal forests [52].
The effect of updating the respiration algorithm on global MR was substantial, increasing the global total by 38%, largely due to a substantial increase in basal respiration. The size of this increase was somewhat surprising in that despite the higher base respiration rate, the updated algorithm includes a mediumterm temperature acclimation response that is not present in the current (default) algorithm, but acclimation had only modest impacts on annual MR. The increase in MR was reflected in large decreases in global NPP (19%) and LAI (15%), due to the smaller fraction of carbon available for allocation. However, these effects led to only a modest reduction in global GPP, which declined by 4%. The updated model tended to lower GPP the most in the tropics, where increases in MR were the highest.
Many of the physiological components in global land surface models are dependent on the latest advances in empirical ecological and plant physiological research. Here, we brought the respiration algorithms of ELM into alignment with recent empirical advances and found significant changes in the global carbon cycle, and its temporal and spatial dynamics.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// zenodo.org/record/5217876.

Funding
This research was supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the US Department of Energy, Office of Science, Office of Biological and Environmental Research (including Grant DE-SC0012677 to P B R), Biological Integration Institutes Grant NSF-DBI-2021898 (to P B R), and by NSF grants IIS-1908104, OAC-1934634, and a grant from C3.ai (to AB).