Land Carbon Budget with Growing Plants & Three Decomposing Pools

CO2 Fertilization:

Relative increase in NPP per doubling of atmopsheric CO2. Experimental value ~ 25%

Nitrogen Fertilization:

Increase in plant carbon when nitrogen demand is completely satisfied. Baseline value ~ 20%

Deforestation & Disturbance:

Deforestation history is prescribed. Adjust the magnitude by scaling the disturbance & recovery fluxes up and down. Baseline value = 2.

Respiration & Decomposition

Plant Growth and Initialization

Model includes deforestation, regrowth, temperature sensitivity to respiration/decomposition, fertilization by CO2 and nitrogen.

Land Carbon Fluxes and Pools

download model output

Assumed timeseries of forcing data

download model input


model diagram

This very simple box model of the terrestrial carbon cycle combines logistic plant growth with decomposing soil organic mater in a series of three cascading pools. The idea is to enhance plant growth rates according to CO2 fertilization and increase the carrying capacity according to nutrient availability, which grows over time due to nitrogen deposition. Plant mortality (conversion to dead litter) is in turn enhanced by deforestation or disturbance. Dead carbon in litter and soil pools is metabolized by microbes whose respiration efficiency is adjustable and depends on temperature.

We track four carbon pools: plants (P), litter (L), “fast” soil (F), and “slow” soil (S). The amount of carbon stored in each pool is determined by four predictive model equations:

\[ \frac{dP}{dt} = NPP - MORT \] \[ \frac{dL}{dt} = MORT - \frac{L}{\tau_L}Q_{10} ^{\frac{T-T_0}{10}} \] \[ \frac{dF}{dt} = \left[ \frac{L}{\tau_L}(1-\epsilon) - \frac{F}{\tau_F} \right] Q_{10} ^{\frac{T-T_0}{10}} \] \[ \frac{dS}{dt} = \left[ \frac{F}{\tau_F}(1-\epsilon) - \frac{S}{\tau_S} \right] Q_{10} ^{\frac{T-T_0}{10}} \]

Symbols used: * NPP = Net primary production (photosynthesis minus plant respiration, in GtC/yr) * MORT = plant mortality (transfers carbon from live plants to dead litter, GtC/yr) * \(\tau_L\), \(\tau_F\), \(\tau_S\) = turnover time for carbon in litter, fast, and slow soil pools (years) * \(\epsilon\) = efficiency of microbial respiration (fraction of what they “eat”" that turns into CO2) * \(Q_{10}\) = fractional increase in decomposition rates per 10 degrees of warming * T = temperature (Celsius), prescribed from a simple climate model using the IPCC SRES A2 emissions scenaro.

Plant growth (NPP) is modeled assuming Logistic growth as widely used in biology and ecology to describe populations growth under environmental constraints. Essentially this is a modification of exponential growth with a growth rate modified (reduced) as the population approaches limits set by nutrients or other environmental conditions.

\[ NPP = g P \left(1 - \frac{P}{NK} \right) \]

It’s helpful to think of \(g(1-\frac{P}{NK})\) as a modified or constrained growth rate, which is the probability of reproduction for each individual in a population per timestep. Here \(N\) is a “nutrient status” which multiplies \(K\), the “carrying capacity” or nutrient-limited equilibrium size of the plant carbon pool. As the plant pool \(P\) approaches the carrying capacity \(NK\), the modified growth rate \(g(1-\frac{P}{NK})\) approaches zero. When \(P\) is very small, the plant pool grows exponentially.

The nutrient status \(N\) multiplies the original (unfertilized) carrying capacity according to a prescribed function that ramps up following an acrtangent curve from 1.0 in the year 1800 to 1.2 in 2150. This factor can be adjusted with a “Nitrogen fertilization” slider in the left-hand control panel of the website.

The growth rate \(g\) of the plants is also enhanced by increasing atmospheric CO2 (CO2 fertilization). This is represented in the model as

\[g = g_0 \left(1 + \beta ln \frac{CO2}{CO2_{init}} \right)\]

where \(g_0\) is the equilibrium growth rate at the initial CO2 concentration \(CO2_{init}\), and \(\beta\) represents the strength of the CO2 fertilization effect on NPP. Analyses of forest inventory data suggest that \(\beta \approx 0.36\) for modern temperate forests. This value corresponds to an increase of NPP by 25% per doubling of atmospheric CO2. The strength of this effect can be adjusted using another model control slider. CO2 is prescribed from a simple climate model using the IPCC SRES A2 emissions scenaro.

Plant mortality \(MORT\) is modeled as a fraction of the living plant carbon pool, plus a prescribed disturbance of the plant pool due to deforestation:

\[ MORT = death {P} + disturbance = \frac{g_0}{\lambda} P + disturbance\]

where \(\lambda\) is the average lifetime of plants in years.

The disturbance rate due to historical deforestation is prescribed in the model using a very simple (1+cosine) function that ramps up from zero in 1825 to a maximum of 2 GtC/yr in 1975, and then back down to zero in 2125. The magnitude of this term can be scaled up or down with yet another model input slider on the website.

The model is initialized in the year 1800 by setting values of the initial plant pool (\(P_{eq}\), default value is 500 GtC) and global NPP (\(NPP_{eq}\), default value is 60 GtC/yr). Carrying capacity in this undisturbed state is given by

\[ K = P_{eq} \left(1 - \frac{1}{\lambda} \right)\].

The equilibrium growth rate of the plant pool for the initial condition is given by

\[ g_0 = \frac{NPP_{eq}}{P_{eq} \left( 1 - \frac{P_{eq}}{K} \right)}\].

The program that does the calculation is very simple, and the code that controls this website is surprisingly simple too! It is all written in the programming language R, using a web programming package called shiny. You can read all about it on the “Website Code” tab to the right.

How this website works (including all the code!)

This website is controlled using the R package “shiny.” There are five important components:

  • Simple program (land.R) that integrates the land carbon cycle model
  • Script ( generates plots of model output (carbon fluxes and pools over time)
  • Script (plot.input.R) generates plots of model inputs (CO2, temperature, disturbance, and nutrient status over time)
  • The user interface (ui.R) that controls the sliders and passes inputs to the model
  • The “server” (server.R) that responds to user interface events in the web browser

Scroll down or click links in the list above to read all about it!

This is the actual box model that calculates land carbon fluxes and pools


land <- function(tau.litter=2.,, tau.slow=500, eff.microbes=0.8, plant.eq=500, NPP.eq=60., betaCO2=0.36, longevity=2, Q10=2.0, n.limitation=0.2, disturb.factor=1){

# Argument list:
# tau.plants : turnover time for live plants (default = 2 yr)
# tau.litter : turnover time for dead plant material (default = 20 yr)
# tau.soil   : turnover time of (armored) soil organic matter (default = 1000 yr)
# eff.microbes : efficiency of microbial respiration 
# fert.factor   : fractional increase in NPP per doubling of CO2

# Initialize a bunch of arrays for later plotting
nYears <- 500
NEE <- replicate(nYears,NA)
NPP <- replicate(nYears,NA) <- replicate(nYears,NA)
resp.litter <- replicate(nYears,NA) <- replicate(nYears,NA)
resp.slow <- replicate(nYears,NA)
plant <- replicate(nYears,NA)
litter <- replicate(nYears,NA)
fast.soil <- replicate(nYears,NA) 
slow.soil <- replicate(nYears,NA)  
mortality <- replicate(nYears,NA)

# Initialize mass of carbon (kg C) in plants, soil, and passive pools
# These are calculated as steady-state solutions to the differential equations
eq.capacity <- plant.eq / (1-1/longevity)  # resource-limited "carrying capacity"
growth.rate  <- NPP.eq / (plant.eq*(1-plant.eq/eq.capacity))
death.rate <- growth.rate/longevity  # fractional death per year

# Initialize each pool (GtC) to be in equilibrium with NPP and decay
plant[1] <- plant.eq 
litter[1] <- tau.litter * death.rate * plant[1] 
fast.soil[1] <- * (1-eff.microbes) * litter[1] 
slow.soil[1] <- tau.slow/ * (1-eff.microbes) * fast.soil[1]

# Read and set up driver data for CO2, temperature, disturbance, & nutrients
# CO2 and temperature from a box model driven with IPCC SRES A2 <- read.table('data/', header=TRUE)
CO2 <-$CO2
temp <-$temp

deforestation <-$deforestation * disturb.factor
abandonment <-$abandonment * disturb.factor
nutrient <-$nutrient * n.limitation

# Start plant capacity at equilibrium value
capacity <- eq.capacity

# Integrate the model
for (i in 1:(nYears-1)){

  # Apply CO2 fertilization to growth (NPP)
  enhanced.growth.rate <- growth.rate * (1.0 + betaCO2 * log(CO2[i]/CO2[1]))

  # Subtract half of deforestation from plant capacity, but add back abandonment
  capacity <- capacity + longevity * (-deforestation[i] + abandonment[i]*2)

  # Enhance plant capacity due to nutrient deposition
  eff.capacity <- capacity * (1 + nutrient[i])

  # Apply nutrient limitation & fertilization to get updated NPP    
  NPP[i] <- enhanced.growth.rate * plant[i] * (1 - plant[i]/eff.capacity)

  # Apply disturbance & mortality                                             
  mortality[i] <- death.rate * plant[i] + deforestation[i]

  # Apply temperature enhancement of respiration
  resp.enhancement <- Q10 ^ ((temp[i]-temp[1])/10)

  # Calculate respiration for each pool          
  resp.litter[i] <- eff.microbes * litter[i]/tau.litter * resp.enhancement[i] <- eff.microbes * fast.soil[i]/ * resp.enhancement
  resp.slow[i] <- slow.soil[i]/tau.slow * resp.enhancement[i] <- resp.litter[i] +[i] + resp.slow[i]

  # Net ecosystem exchnage is repiration minus NPP
  NEE[i] <-[i] - NPP[i]

  # Update all the carbon pools in plants and soils
  plant[i+1] <- plant[i] + NPP[i] - mortality[i]
  litter[i+1] <- litter[i] + mortality[i] - litter[i]/tau.litter * resp.enhancement
  fast.soil[i+1] <- fast.soil[i] + resp.enhancement * (
                    (1.-eff.microbes) * litter[i]/tau.litter -
  slow.soil[i+1] <- slow.soil[i] + resp.enhancement * (
                    (1.-eff.microbes) * fast.soil[i]/ -


return(list(parameters=list(tau.litter=tau.litter,, tau.slow=tau.slow, 
                            eff.microbes=eff.microbes, plant.eq=plant.eq, NPP.eq=NPP.eq, 
                            betaCO2=betaCO2, longevity=longevity, Q10=Q10,
                            n.limitation=n.limitation, disturb.factor=disturb.factor),
           output=data.frame(NEE=NEE, plant=plant, litter=litter, fast.soil=fast.soil,
                      slow.soil=slow.soil, NPP=NPP, CO2=CO2, deforestation=deforestation,
                      abandonment=abandonment, resp.litter=resp.litter,, 
                      resp.slow=resp.slow,, temp=temp, 
                      nutrient=nutrient, mortality=mortality)


This is the script that plots model output fluxes and pools <- function(history){

  years <- 1800:2299

  orig.par <- par(no.readonly=TRUE)
  par(mfrow=c(3,1), cex.lab=1.5, cex.axis=1.5, cex.main=2)

  # Top panel is NEE
  plot(years, NEE, main='Net Ecosystem Flux to Atmosphere', ylab='GtC/yr', 
       type='l', lwd=6)

  # Second panel is Gross fluxes
  plot(years, NPP, main='Gross Fluxes', ylab='GtC/yr', col='darkgreen', 
       type='l', lwd=6, ylim=plot.limits)
  lines(years,, col='red', type='l', lwd=6)
  legend('right', c('NPP, RESP'), col=c('darkgreen','red'), lwd=c(6,6))

  # Bottom panel is carbon pools
  plot.limits <- c(0,max(c(plant, litter, fast.soil, slow.soil),na.rm=TRUE))
  plot.col = c('darkgreen','darkolivegreen2','burlywood','chocolate4')
  plot(years, plant, main='Carbon Pools', ylab='GtC', col=plot.col[1], 
  lines(years, litter, col=plot.col[2], type='l', lwd=6)
  lines(years, fast.soil, col=plot.col[3], type='l', lwd=6)
  lines(years, slow.soil, col=plot.col[4], type='l', lwd=6)
  legend('topright', c('Plants','Litter','Fast','Slow'), col=plot.col, 


This is the script that plots model output inputs


plot.input <- function(history){

  years <- 1800:2299

  orig.par <- par(no.readonly=TRUE)
  par(mfrow=c(4,1), cex.lab=1.5, cex.axis=1.5, cex.main=2)

  # Top panel is CO2
  plot(years, history$CO2, main='Atmospheric CO2', ylab='ppmv', 
       type='l', lwd=6)

  # Second panel is nutrient status
  plot(years, history$nutrient, main='Nutrient Status', ylab='ratio', col='darkgreen', 
       type='l', lwd=6)

  # Third panel is deforestation & reforestation
  plot(years, history$deforestation, main='Disturbance & Recovery', ylab='GtC/yr', 
       col='red', type='l', lwd=6)
  lines(years, history$abandonment, col='green',lwd=6)
  legend('topright', c('Deforestation','Reforestation'), col=c('red','green'),
         lwd=c(6,6), cex=0.8)

  # Fourth panel is temperature 
  plot(years, history$temp, main='Warming', ylab='Celsius', col='red', 
       type='l', lwd=6)


This is the user interface that actually controls this website



# Define user interface for land carbon model

  #  Application title
  headerPanel("Land Carbon Budget with Growing Plants & Three Decomposing Pools"),

  # Sidebar on the left with controls for the user to specify atmospheric properties

    # Adjust strength of CO2 fertilization
      img(src='CO2.png', align='left'),
      h4('CO2 Fertilization: '),
      p('Relative increase in NPP per doubling of atmopsheric CO2.
        Experimental value ~ 25%'),
      sliderInput("betaCO2", "(percent)", 
                  min = 1, max = 100, value = 25, step=1)

    # Adjust strength of nitrogen fertilization
      img(src='fertilizer.png', align='left'),
      h4('Nitrogen Fertilization: '),
      p('Increase in potential NPP when nitrogen demand 
        is completely satisfied. Baseline value ~ 20%'),
      sliderInput("n.limitation", "(percent:)", 
                  min = 1, max = 50, value = 20, step=1)

    # Adjust strength of deforestation/disturbance
      img(src='deforestation.png', align='left'),
      h4('Deforestation & Disturbance:'),
      p('Deforestation history is prescribed. Adjust the magnitude by scaling  
        the disturbance & recovery fluxes up and down. Baseline value = 2.'),
      sliderInput("disturb.factor", 'GtC/yr in 1975',
                  min = 0.1, max = 3, value = 2, step= .1)

    # Adjust respiration/decomposition parameters
      h4('Respiration & Decomposition'),
      checkboxInput('adjustResp', 'Adjust Parameters'),
        sliderInput("Q10", "Temperature sensitivity of respiration (Q10).
                    Experimental value ~ 2", 
                    min = 1, max = 3, value = 2, step=.1),
        sliderInput("eff.microbes", "Efficiency of Microbial Respiration (percent)", 
                    min = 10, max = 95, value = 80, step=1),
        sliderInput("tau.litter", "Litter turnover time (years)", 
                    min = 1, max = 20, value = 2, step= 1),
        sliderInput("", "Fast Soil turnover time (years)", 
                    min = 1, max = 50, value = 5, step= 1),
        sliderInput("tau.slow", "Slow turnover time (years)", 
                    min = 100, max = 2000, value = 600, step=10)   

    # Adjust plant growth parameters
      h4('Plant Growth and Initialization'),
      checkboxInput('adjustPlants', 'Adjust Parameters'),
        sliderInput("longevity", "Average Lifetime of Plants (years)", 
                    min = 1, max = 20, value = 2, step=1),
        sliderInput("plant.eq", "Baseline Plant Carbon Pool (GtC)", 
                    min = 100, max = 1000, value = 500, step=10),
        sliderInput("NPP.eq", "Baseline NPP (GtC/yr)", 
                    min = 10, max = 100, value = 60, step=1)

  # Show a tabset in the main panel of the browser that displays model output
    h4('Model includes deforestation, regrowth, temperature sensitivity
       to respiration/decomposition, fertilization by CO2 and nitrogen.'),


      # First tab shows plots of model output for both control and modified atmosphere
      tabPanel("Model Output", 
               h4("Land Carbon Fluxes and Pools"),
               tags$style(type="text/css", ".tab-content { overflow: visible; }")), 

      # Second tab shows plots of model output for both control and modified atmosphere
      tabPanel("Model Input", 
               h4("Assumed timeseries of forcing data"),
               tags$style(type="text/css", ".tab-content { overflow: visible; }")), 

      # Third tab displays a brief model description
      tabPanel("Model Description", 

      # Fourth tab displays website code
      tabPanel("Website code", 

This “server” script responds to user interface events and creates browseer content



# Source required R scripts

shinyServer(function(input, output) {

  # Run the model, saving output to a structured list called "modelOutput"
  modelOutput <- reactive({
         tau.slow=input$tau.slow, eff.microbes=input$eff.microbes/100, 
         plant.eq=input$plant.eq, NPP.eq=input$NPP.eq, 
         longevity=input$longevity, betaCO2=input$betaCO2/log(2)/100., 
         Q10=input$Q10, n.limitation=input$n.limitation/20, 

  # Ladder plot of carbon fluxes and pools
  output$land.plot <- renderPlot({
  }, height=800)

  # Ladder plot of model inputs
  output$input.plot <- renderPlot({
  }, height=800)