Skip to contents


Hector’s land system is divided into three pools: vegetation (\(C_v\)), detritus (\(C_d\)), and soil (\(C_s\)); a fourth pool, earth (\(C_earth\)) is used only for anthropogenic emissions and uptake. These pools are global by default, but they can be divided into \(n\) “biomes” with unique parameters and initial pool sizes.

The total C flux from atmosphere to land at time \(t\) (\(F_{L}(t)\); also known as Net Biome Production, NBP) is the net of C sequestration by photosynthesis (net primary productivity, \(NPP_i(t)\)) and C release by heterotrophic respiration (\(RH_i(t)\)) summed across each biome \(i\), and global land use change (\(LUC(t)\)) (note that this calculation does not currently include fire, another important land C source):

\[ F_{L}(t) = \sum_{i=1}^{n} \left( NPP_i(t) - RH_i(t) \right) - LUC(t) \]

NPP is calculated as the product of a user-specified pre-industrial value (\(NPP_0\), default is 56.2 Pg C/year), a CO2 fertilization multiplier (\(f(C_{atm}, \beta)\)), and the effect of land use change (LUC) on vegetation carbon stocks:

\[ NPP_i(t) = NPP_{0,i} \times f(C_{atm}, \beta_i) \times f(LUC_{v}) \]

The CO2 fertilization multiplier is a function of the current atmospheric CO2 concentration (\(C_{atm}\)), the initial (“pre-industrial”) CO2 concentration (\(C_0\)), and biome-specific CO2 fertilization parameter \(\beta_i\) (default = 0.36):

\[ f(C_{atm}, \beta_i) = 1 + \beta_i \left( \log(\frac{C_{atm}}{C_0})\right) \]

The LUC multiplier accounts for the effect of vegetation losses or gains on NPP (although we recognize that LUC does not always result in ‘lost’ NPP); this behavior is new as of Hector v3. It is computed as a fraction based on preindustrial (end-of-spinup) vegetation carbon (\(C_{v}\)) and the running sum of LUC-driven losses from – or gains to – vegetation:

\[ f(LUC_{v}) = \frac{C_{0v} - \sum LUC_{veg}}{C_{0v}} \]

Heterotrophic respiration occurs in soil (\(RH_s(t)\)) and detritus (\(RH_d(t))\)) as a function of biome-specific atmospheric temperature anomaly (\(T_i\), calculated as the global atmospheric temperature anomaly scaled by a biome-specific warming factor \(\delta\)) and the size of the current soil and detritus C pools (\(C_s\) and \(C_d\), respectively).

\[ T_i(t) = \delta \, T_G(t) \] \[ RH_s = \frac{1}{50} C_s Q_{10}^{T_{i}(t) / 10}\] \[ RH_d = \frac{1}{4} C_d Q_{10}^{T_{i}(t) / 10}\]

Note that the residence times of soil and detritus C pools – 50 and 4 years, respectively – are currently hard-coded values in Hector.

Fluxes associated with land use change (\(LUC\)) are prescribed via the luc_emissions variable in the INI file. Example time series associated with shared socioeconomic pathways (SSPs) are included with Hector in the corresponding inst/input/tables/ssp**_emiss-constraints_rf.csv files. There are two LUC fluxes, land uptake and loss, and both add or remove carbon to the land pools in proportion to those pools’ size. (This behavior is new as of version 3 of the model.) Note that the cumulative \(LUC_{v}\), LUC effects on vegetation, affects NPP as documented above.

In each biome at each time step, the vegetation pools change as follows (biome, \(i\) and time, \(t\) subscripts omitted for clarity): C gain from NPP is distributed to vegetation, detritus, and soil according to fraction parameters \(f_{nv}\), \(f_{nd}\), and \(f_{ns}\), respectively. Similarly, each pool bears a fraction of the total C loss from land-use change (\(F_{LC}\)), again according to fractions \(f_{lv}\), \(f_{ld}\), and \(f_{ls}\). The detritus and soil pools also lose C via heterotrophic respiration (\(RH_d\) and \(RH_s\), respectively). Finally, a fixed fraction of vegetation C is transferred to detritus as litterfall (\(f_{vd}\)), and a fraction of both vegetation (\(f_{vs}\)) and detritus (\(f_{ds}\)) C are transferred to soil.

\[ \frac{dC_v}{dt} = NPP\, f_{nv} - C_v (f_{vd} + f_{vs}) - F_{LC} f_{lv} \] \[ \frac{dC_d}{dt} = NPP\, f_{nd} + C_v f_{vd} - C_d f_{ds} - RH_d - F_{LC} f_{ld} \] \[ \frac{dC_s}{dt} = NPP\, f_{ns} + C_v f_{vs} - C_d f_{ds} - RH_s - F_{LC} f_{ls} \]

Anthropogenic changes to the earth system are driven by fossil fuel and industrial emissions (FFI) as well as, in some scenarios, direct air carbon capture and storage (DACCS). These fluxes are occur out of, or into, an ‘earth’ pool that is used for mass-balance checks:

\[ \frac{dC_{earth}}{dt} = DACCS - FFI\]

Finally, the atmosphere changes as a result of all the fluxes into or out of it: anthropogenic emissions and uptake; land use change emissions and uptake; fluxes into or out of the ocean; and heterotrophic respiration and net primary production.

\[ \frac{dC_{atm}}{dt} = FFI - DACCS + LUC_{e} - LUC_{u} + O_{i} - O_{u} + RH - NPP \]

The entire terrestrial C cycle is summarized in the following diagram:

Summary of terrestrial C cycle in Hector. The full Hector C cycle also includes ocean and “Earth” C pools that exchange C with the atmosphere, but these are omitted for clarity.


The state of the global C cycle in Hector is defined by a length 6 vector c, which describes the amount of C (in Pg C) in the atmosphere (0), vegetation (1), detritus (2), soil (3), ocean (4), and “Earth” (5) at the current time step. At each Hector time step, c is solved by the CarbonCycleSolver::run() function (defined in src/carbon-cycle-solver.cpp) in the following general steps:

  1. Retrieve the current estimates of various pools from variables defined by the C cycle model and assign them to the relevant parts of c (getCValues).
  2. Integrate from the beginning of the time step using updated slow parameters (slowparameval).
  3. Solve the ordinary differential equations associated with the C cycle model based on its time derivatives (defined in calcderivs).
  4. Re-calculate C cycle model internal variables based on the solved values of c (stashCValues).
  5. Record the model state (record_state).

These steps are described in more detail below. The complete code for the terrestrial C cycle lives in src/simpleNbox.cpp file (with associated headers in inst/include/simpleNbox.hpp).

Setting the c vector (getCValues)

This function sets the values of the c vector to the corresponding state variables in the C cycle model.

  • The atmospheric C pool for the current time step is defined by atmos_co2, which is always a scalar unit quantity (because CO2 is assumed to be well-mixed in the atmosphere).

  • Vegetation, detritus, and soil pools at the current time step are stored in veg_c, detritus_c, and soil_c, respectively. These are vector quantities (“stringmaps”), with one value for each biome (the default biome name is “global”), so the value stored in the corresponding slot in c is the sum across all biomes. The ocean C pool is retrieved by a call to the current ocean model’s getCValues function. (NOTE: This means that simpleNBox model depends on the existence of an ocean C cycle model.)

  • Finally, the “Earth” C pool is stored in earth_c, which is also a scalar unit quantity.

All of these variables have are in units Pg C.

Slow parameter evaluation (slowparameval)

  • Perform the ocean C model’s slow parameter evaluation (omodel->slowparameval)

  • Calculate CO2 fertilization effect for each biome (co2fert[biome]) using the function calc_co2fert

    • …unless you are in the spin-up phase, in which case assume no effect (co2fert = 1)
  • Calculate temperature effect on respiration of detritus (tempfertd) and soil (tempferts).

    • Note that soil warming is “sticky” – it can only increase, not decline (even if the temperature declines between time steps)
    • Note also that while detritus warming is based on the actual current biome temperature, the soil warming is based on a 200-year running mean

Note that these effect parameters are calculated but not yet applied. They will be applied to the actual NPP and RH estimates during the ODE solving loop.

C cycle derivatives (calcderivs)

  • Calculate NPP for each biome (npp function)
  • Partition NPP C gain into vegetation, detritus, and soil pools according to fractions
  • Calculate RH for detritus and soil
  • Partition C losses from RH according to fractions
  • Calculate and partition litter flux
  • Store changes in pools in dcdt vector (same structure as c)

Storing the results (stashCValues)

  • Reset the atmospheric CO2 state variable (atmos_co2) to the current value from c
  • Calculate the difference between vegetation, detritus, and soil C calculated by the ODE solver (i.e. in c) and what was stored in the current state variables (veg_c, detritus_c, soil_c). Store these differences in veg_delta, det_delta, soil_delta.
  • Apply these additional fluxes to each biome, weighing by the biome’s current gross primary productivity (NPP + RH). I.e.:
    • Calculate the total GPP across all biomes (npp_rh_total = sum_npp() + sum_rh())
    • Calculate each biome’s weight (wt[biome] = (npp(biome) + rh(biome)) / npp_rh_total)
    • Add the corresponding difference multiplied by the weight to the biome’s current pool (e.g. veg_c[biome] = veg_c[biome] + veg_delta * wt[biome])
    • NOTE: This is a workaround. In the future, Hector will use the ODE solver to separately solve all of the boxes in the multi-biome system.
  • Stash the ocean C pools using the omodel->stashCValues function
  • Stash the Earth C pool (earth_c)
  • Check for mass conservation by take the sum of c (total C in system) and comparing it against the sum from the previous timestep (masstot).
    • NOTE: This check is skipped the first time this function is evaluated, at which point masstot is initialized to the calculated sum.
  • If Hector is running under atmospheric CO2 constraint, re-calculate the atmospheric CO2 pool to match the constraint and move any residual C to/from the deep ocean pool.

Recording the state (record_state)

This just sets the time-series versions of all of the C cycle state variables (e.g. atmos_c_ts, veg_c_tv, detritus_c_tv, tempfertd_tv) to their corresponding state variables at the current time step (e.g. atmos_c, veg_c, detritus_c, tempfertd).

Setting (setData) and retrieving (getData) values

The following values related to the terrestrial C cycle can be set via setData.

Description Unit C++ macro R function C++ variable Needs date Biome-specific Set Get
Initial atmospheric CO2 pool PgC D_ATMOSPHERIC_CO2 ATMOSPHERIC_CO2() C0 No No Yes Yes
Pre-industrial CO2 concentration ppmv CO2 D_PREINDUSTRIAL_C02 PREINDUSTRIAL_CO2() C0 No No Yes Yes
Atmospheric C constraint residual ppmv CO2 D_ATMOSPHERIC_C_RESIDUAL residual No No No Yes
Vegetation C pool PgC D_VEGC veg_c No Yes Yes Yes
Detritus C pool PgC D_DETRITUSC detritus_c Yes Yes Yes Yes
Soil C pool PgC D_SOILC soil_c No Yes Yes Yes
Radiative forcing from albedo unitless D_RF_ALBEDO RF_ALBEDO() RF_albedo Yes No Yes Yes
Frac. NPP to vegetation unitless D_F_NPPV F_NPPV() f_nppv No Yes Yes Yes
Frac. NPP to detritus unitless D_F_NPPD F_NPPD() f_nppd No Yes Yes Yes
Frac. veg. C to litter unitless D_F_LITTERD F_LITTERD() f_litterd No Yes Yes Yes
Initial NPP flux PgC / year D_NPP_FLUX0 npp_flux0 No Yes Yes Yes
Fossil fuel CO2 emissions PgC / year D_FFI_EMISSIONS FFI_EMISSIONS() ffiEmissions Yes No Yes Yes
Direct air carbon capture and storage PgC / year D_DACCS_UPTAKE DACCS_UPTAKE() daccsUptake Yes No Yes Yes
LUC emissions PgC / year D_LUC_EMISSIONS LUC_EMISSIONS() lucEmissions Yes No Yes Yes
Atmospheric CO2 constraint ppmv CO2 D_CO2_CONSTRAIN CO2_constrain Yes No Yes No
CO2 fertilization factor unitless D_BETA BETA() beta No Yes Yes Yes
Biome warming multiplier unitless D_WARMINGFACTOR warmingfactor No Yes Yes Yes
RH temperature sensitivity unitless D_Q10_RH Q10_RH() q10_rh No Yes Yes Yes

Other routines

  • init() – Initialize the C cycle model
    • Set initial values for the default biome ("global") for the following variables:
      • co2fert (CO2 fertilization multiplier) – 1.0
      • warmingfactor (biome warming multiplier relative to global average) – 1.0
      • residual (difference between modeled atmospheric CO2 and constraint) – 0
      • tempfertd, tempferts (Q10 warming effects on detritus and soil respiration, respectively) – 1.0
    • Add "global" to the biome list
    • Register C cycle model capabilities, dependencies, and inputs
  • reset(time) – Resets the state of the C model to time time.
    • Set the values of all current-time state variables (e.g. atmos_co2, veg_c, tempfertd) to their time-series values at time time (e.g. atmos_c_ts.get(time), veg_c_tv.get(time), tempfertd_tv.get(time))
    • Re-calculate the CO2 fertilization effect for each biome (co2fert). If in spinup, force this to 1.0.
    • “Truncate” all time series variables up to the current time; i.e. erase all of the values after time
    • Set tcurrent to time
  • sanitychecks() – Ensures the validity of the current model state. Specific checks are:
    • All pools and fluxes (for each biome) must be greater than or equal to 0
    • Partitioning fraction parameters (e.g. f_nppv, f_litterd) must be between 0 and 1, and must sum to 1
  • prepareToRun() – Perform final checks and variable initialization before running a simulation
    • Check that all biome-specific pools, fluxes, and parameters have data for all biomes. (The only optional biome-specific parameter is warmingfactor, which is set to 1.0 if missing). Also check that biome-specific CO2 fertilization (beta) and Q10 (q10_rh) parameters are valid.
    • If no albedo forcing data are provided, set the value to the MAGICC default of -0.2
    • Initialize the atmospheric CO2 concentration (CO2_concentrations, ppm) and C pool (atmos_co2, PgC) to the user-specified initial value (C0; converted accordingly to PgC for CO2_concentrations)
    • Run the checks in sanitychecks()