Terrestrial C cycle

Science

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 CO₂ 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 CO₂ fertilization multiplier is a function of the current atmospheric CO₂ concentration (\(C_{atm}\)), the initial (“pre-industrial”) CO₂ concentration (\(C_0\)), and biome-specific CO₂ 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.

Implementation

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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ constraint, re-calculate the atmospheric CO₂ 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 (CO₂ fertilization multiplier) – 1.0
    • warmingfactor (biome warming multiplier relative to global average) – 1.0
    • residual (difference between modeled atmospheric CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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()