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.

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 (plot.land.R) 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!

**land.R:**

land <- function(tau.litter=2., tau.fast=20., 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)
resp.total <- replicate(nYears,NA)
resp.litter <- replicate(nYears,NA)
resp.fast <- 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] <- tau.fast/tau.litter * (1-eff.microbes) * litter[1]
slow.soil[1] <- tau.slow/tau.fast * (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
driver.data <- read.table('data/driver.data.txt', header=TRUE)
CO2 <- driver.data$CO2
temp <- driver.data$temp
deforestation <- driver.data$deforestation * disturb.factor
abandonment <- driver.data$abandonment * disturb.factor
nutrient <- driver.data$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
resp.fast[i] <- eff.microbes * fast.soil[i]/tau.fast * resp.enhancement
resp.slow[i] <- slow.soil[i]/tau.slow * resp.enhancement
resp.total[i] <- resp.litter[i] + resp.fast[i] + resp.slow[i]
# Net ecosystem exchnage is repiration minus NPP
NEE[i] <- resp.total[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 -
fast.soil[i]/tau.fast)
slow.soil[i+1] <- slow.soil[i] + resp.enhancement * (
(1.-eff.microbes) * fast.soil[i]/tau.fast -
slow.soil[i]/tau.slow)
}
return(list(parameters=list(tau.litter=tau.litter, tau.fast=tau.fast, 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.fast=resp.fast,
resp.slow=resp.slow, resp.total=resp.total, temp=temp,
nutrient=nutrient, mortality=mortality)
))
```

}

**plot.land.R:**

```
plot.land <- function(history){
attach(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)
grid(col='black')
# Second panel is Gross fluxes
plot.limits=range(c(NPP,resp.total),na.rm=TRUE)
plot(years, NPP, main='Gross Fluxes', ylab='GtC/yr', col='darkgreen',
type='l', lwd=6, ylim=plot.limits)
lines(years, resp.total, col='red', type='l', lwd=6)
legend('right', c('NPP, RESP'), col=c('darkgreen','red'), lwd=c(6,6))
grid(col='black')
# 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],
ylim=plot.limits)
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,
lwd=rep(6,4))
grid(col='black')
detach(history)
par(orig.par)
}
```

**plot.input.R:**

```
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)
grid(col='black')
# Second panel is nutrient status
plot(years, history$nutrient, main='Nutrient Status', ylab='ratio', col='darkgreen',
type='l', lwd=6)
grid(col='black')
# 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)
grid(col='black')
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)
grid(col='black')
par(orig.par)
}
```

**ui.R:**

```
library(shiny)
library(markdown)
source('model/land.R')
source('model/plot.land.R')
source('model/plot.input.R')
# Define user interface for land carbon model
shinyUI(pageWithSidebar(
# 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
sidebarPanel(
# Adjust strength of CO2 fertilization
wellPanel(
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
wellPanel(
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
wellPanel(
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
wellPanel(
h4('Respiration & Decomposition'),
checkboxInput('adjustResp', 'Adjust Parameters'),
conditionalPanel(condition="input.adjustResp",
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("tau.fast", "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
wellPanel(
h4('Plant Growth and Initialization'),
checkboxInput('adjustPlants', 'Adjust Parameters'),
conditionalPanel(condition="input.adjustPlants",
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
mainPanel(
h4('Model includes deforestation, regrowth, temperature sensitivity
to respiration/decomposition, fertilization by CO2 and nitrogen.'),
tabsetPanel(
# First tab shows plots of model output for both control and modified atmosphere
tabPanel("Model Output",
h4("Land Carbon Fluxes and Pools"),
plotOutput('land.plot'),
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"),
plotOutput('input.plot'),
tags$style(type="text/css", ".tab-content { overflow: visible; }")),
# Third tab displays a brief model description
tabPanel("Model Description",
includeHTML('doc/model.description.html')),
# Fourth tab displays website code
tabPanel("Website code",
includeMarkdown('doc/website.code.md'))
)
)
))
```

**server.R:**

```
library(shiny)
# Source required R scripts
source('model/land.R')
source('model/plot.land.R')
shinyServer(function(input, output) {
#
# Run the model, saving output to a structured list called "modelOutput"
modelOutput <- reactive({
land(tau.litter=input$tau.litter, tau.fast=input$tau.fast,
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,
disturb.factor=input$disturb.factor/2)
})
# Ladder plot of carbon fluxes and pools
output$land.plot <- renderPlot({
plot.land(modelOutput())
}, height=800)
# Ladder plot of model inputs
output$input.plot <- renderPlot({
plot.input(modelOutput())
}, height=800)
})
```