One may assume that the general dynamics of the coupled C and N cycling in terrestrial ecosystems can be described by a small list of processes (functions). The simplest CN-model (“sCN-model”) embodies the following relationships:

  • The plant (divided into an above and a belowground compartment) balances the acquisition of C and N in order to satisfy its stoichiometric requirements (fixed: \(r_{C:N}\)), imposed by new growth (balanced growth condition). It tends to achieve a balance between supply of N through acquisition pathways and the demand imposed by new growth and its respective C:N ratio. \[ {N_{\text{acq}}}(C_r) = y/{r_{\text{C:N}}}[ P(C_l) - r(C_l+C_r) ] \; \; \; \; \; \; (1) \] \(P\) is the C assimilation rate and is a function of leaf mass (\(C_l\)), \(y\) is the growth efficiency, \(r\) a respiration coefficient, implying that plant respiration scales with its size, \(C_l\) and \(C_r\) are the leaf carbon and root pool sizes respectively, and \(N_{\text{acq}}\) is the the N acquisition flux and is a function of the root mass (\(C_r\)).

  • The size of the above-ground pool \(C_l\) (each pool consists of C and N) determines the assimilation rate of C (\(P\)) and has declining marginal returns towards increasing \(C_l\). \[ P(C_l) = I \varepsilon \left(1 - e^{-k_b \sigma C_{\text{l}}} \right) \] Here, \(I\) is PPFD, \(\varepsilon\) is the light use efficiency, \(k_b\) is the light extinction parameter, and \(\sigma\) is the specific leaf area. To simplify the mathematics, I’ve tried to model this with Michaelis-Menten as an alternative: \[ P(C_l) = I \varepsilon \frac{C_l}{C_l + K_P} \]

  • The size of the below-ground pool (\(C_r\)) determines the fraction of net N mineralisation \({N_{\text{min}}}\) that is acquired by the plant (\({N_{\text{acq}}}\)). As for \(P\), the function \(f(C_r)\) represents declining marginal returns in \({N_{\text{acq}}}\) with increasing \(C_r\): \[ {N_{\text{acq}}}= f(C_r) {N_{\text{min}}}\\ f(C_r) = (1 - u) (1 - e^{-k_r C_r } ) \] Here, \(u\) is the unaccessible fraction of net mineralisation, i.e. N losses that are not subject to plant control. Again, in order to simplify mathematics, I’ve tried to model this with Michaelis-Menten as an alternative: \[ f(C_r) = (1-u) \frac{C_r}{C_r + K_N} \]

Example

For given \(I\), \(\varepsilon\), \(N_{\text{min}}\), a balance of above-to-belowground pool size \(a = C_l:C_r\) can be found that satisfies Eq. (1). The imbalance can be expressed as \[ f((1-a)C){N_{\text{min}}}= y/{r_{\text{C:N}}}\; [ P(aC) - rC ] \] Here, \(C=C_l+C_r\) and \(a=C_l/C\). Below is an example plot for the imbalance term as a function of \(a\).

Steady-state solution

In steady state, the total stock of N in the soil is constant. This implies that net mineralisation equals total N inputs through litterfall and N deposition (assuming deposited N bypasses plants). For the sake of simplicity, N fixation and resorption are ignored here. \[ {N_{\text{min}}}= \frac{C}{{r_{\text{C:N}}}\tau} + {N_{\text{in}}}\] Here, \(\tau\) is the residence time of C and N in the plant pool; \({r_{\text{C:N}}}\) is the C:N ratio of plant biomass; and \({N_{\text{in}}}\) is the external input of N, subsuming atmospheric deposition and weathering. Using this, the steady-state solution of the coupled soil-plant system can be found, without relying on prescribed net mineralisation rates, with the balanced-growth condition from above: \[ f((1-a)C)(\frac{C}{{r_{\text{C:N}}}\tau} +{N_{\text{in}}}) = y / {r_{\text{C:N}}}\; [ P(aC) - rC ]\; \; \; \;\; \; \; \;(2) \] and the plant C blance (C assimilation equals respiration plus new growth, replacing litterfall in steady-state): \[ P(aC) - (\tau/y + r)\;C = 0 \; \; \; \; \; \; \; \;(3) \] This set of two equations should in theory allow us to cancel \(C\) and solve for \(a\) as a function of environmental conditions: \(a = f(\varepsilon, I, {N_{\text{in}}})\). Note that \(\varepsilon\) is not treated here as a parameter because it’s itself a function of environmental conditions. Let’s not try to get the analytical solution here. It will involve Weibull functions anyways, and these have to be solved numerically. Let’s just throw it at nleqslv().

library(nleqslv)

## Equation system as a function returning a vector of two values, both of which are to be set to zero by the optimiser
fn <- function( par ){
  eff  <- 0.7
  resp <- 0.1
  r_cton_plant <- 30
  tau_plant <- 10
  n_in <- 1.0

  out <- numeric(2)
  ctot <- par[1]
  f_ag <- par[2]
  
  out[1] <- eval_cnbalance( f_ag, ctot, ( ctot/(r_cton_plant * tau_plant) + n_in), lue, ppfd )
  out[2] <- prod( f_ag * ctot, ppfd, lue ) - ctot * ( tau_plant/eff + resp )
  
  return(out)
}

steadystate <- nleqslv( c(500, 0.5), fn, control = list( allowSingular=TRUE ) )

## test
ctot <- steadystate$x[[1]]
f_ag <- steadystate$x[[2]]
  eff  <- 0.7
  resp <- 0.1
  r_cton_plant <- 30
  tau_plant <- 10
  n_in <- 1.0
  print(paste("should be zero:", eval_cnbalance( f_ag, ctot, ( ctot/(r_cton_plant * tau_plant) + n_in), lue, ppfd ) ))
## [1] "should be zero: 1.84608450057872e-14"
  print(paste("should be zero:", prod( f_ag * ctot, ppfd, lue ) - ctot * ( tau_plant/eff + resp ) ))
## [1] "should be zero: -2.19709934079091e-11"

My problem

Something is apparently not right here. I’ve tried to solve this system of two non-linear equations numerically but the solution I get is nonsense. A problem might be that I cannot constrain the solutions to a (phyically meaningful) interval, \(0<a<1\) in this case. But the nature of the problem shouldn’t require this anyways. So I am asking myself if the second equation above (Eq. 3) is appropriate or if something else is wrong.

So far, I haven’t attempted to analytically solve this. Feel free if you like…

Dynamic CN-model

Above equations can also be implemented in a simple dynamical model that additionally accounts for soil C and N turnover. The results from a spinup to equilibrium look like this (suggesting that there is a physically meaningful solution). In this case (last plot at the very bottom), the ratio of \(C_l/C\) is around 0.83.

## [1] "Steady state f_ag (fraction of aboveground to total plant C):"
## [1] NA
## [1] "Steady state total plant C:"
## [1] NA
## [1] "Steady state net mineralisation:"
## [1] 1.000001
## [1] "Expected steady state net mineralisation:"
## [1] 1.546038