Here the state is $G$. The rate equation can be re-written as “something times $G$”:

$$\frac{dG}{dt} = gc_{in} * F - G * \frac{F}{V}$$

When the inflow $gc_{in}$ equals 0, the first part can be ignored. Then it is easy to see that the time coefficient for $G$, $TC_G$, equals $\frac{V}{F}$, or 1 divided by “what multiplies with the state”. This also looks like an negative exponential equation!

It can also be shown that this solution is correct with the general equation for the time coefficient and the equilibrium value for the state. The equilibrium value for $G$:

$$\frac{dG}{dt}=0 \\ G_{eq} * \frac{F}{V} = gc_{in} * F \\ G_{eq} = gc_{in} * F * \frac{V}{F} \\$$

When using this value in the equation for the time coefficient:

$$TC_G = \huge | \small \frac{gc_{in} * V - G}{gc_{in} * F - G * \frac{F}{V}} \huge | \small = \huge | \small \frac{gc_{in} * V - G}{\frac{F}{V} * (gc_{in} * V - G)} \huge | \small = \huge | \small \frac{1}{\frac{F}{V}} \huge | \small = \frac{V}{F} \\ TC_G = \frac{V}{F} = \frac{4532 × 10^9 \text{ m}^3}{265 \text{ m}^3 \text{ s}^{-1}} = 17.1 × 10^9 \text{ seconds} = 542.30 \text{ years}$$

Note that this TC value is based on the assumption that the lake is perfectly mixed and that $G$ is only removed by outflow (so no decomposition or net precipitation on the lake floor).

This is an negative exponential system, most easy to see when inflow equals 0.

A negative feedback means stability i.e. the system returns to the same equilibrium when a small change in the value of the state is made. A positive feedback is instable, and a small change will move the system away from an equilibrium. This system has a negative feedback. When the amount of $G$ is changed, it returns to the same equilibrium $G$, equal to $gc_{in} * V$. This because the outflow rate increases when $G$ increases, so reducing the amounts in the lake, and the outflow rate decreases when $G$ decreased, increasing the amounts in the lake when the inflow remains unchanged.

The implementation of CUVET1:

# Define a function that returns a named vector with initial state values
InitCUVET1 <- function(W = 0){
  # Gather function arguments in a named vector
  y <- c(W = W)
  
  # Return the named vector
  return(y)
}

# Define a function that returns a named vector with parameter values
ParmsCUVET1 <- function(Volume = 0.0001, # volume of the chamber, m3
                        Area   = 0.01,   # area of leaves, m2
                        Vm   = 0.001,    # Maximum rate of photosynthesis VM in g of CO2 m-2 s-1
                        K      = 1.0,    # K is the M-M constant of photosynthesis
                                         # with respect to CO2 in same units as C
                        Cini   = 10.0){  # initial CO2 concentration in g CO2 m-3
  # Gather function arguments in a named vector
  y <- c(Volume = Volume, 
         Area = Area, 
         Vm = Vm, 
         K = K, 
         Cini = Cini)
  
  # Return named vector with values
  return(y)
}

# Define the function returning the rates of change of the state variables with respect to time
RatesCUVET1 <- function(t, y, parms){
  with(as.list(c(y, parms)),{
    # Compute auxiliary/helper functions

    # Concentration of CO2 inside the chamber
    C   <- W / Volume  
    # Net rate of photosynthesis (no respiration assumed)
    P <- Vm * C / (K + C)
    
    # Compute rate equations
    
    # Net rate of change of amount of CO2 in chamber
    dWdt <- - P * Area
    
    # Gather all rates of change in a vector
    # - the rates should be in the same order as the states (as specified in 'y')
    # - it can be a named vector, but does not need to be
    RATES <- c(dWdt = dWdt)
    
    # Optional: get in/out flow used to compute mass balances (or set MB <- NULL)
    # not included here (thus here use MB <- NULL), see template 3
    MB <- NULL

    # Optional: gather auxiliary variables that should be returned (or set AUX <- NULL)
    # - this should be a named vector or list!
    AUX <- c(C = C, P = P)
    
    # Return result as a list
    # - the first element is a vector with the rates of change (in the same order as 'y')
    # - all other elements are (optional) extra output, which should be named
    outList <- list(c(RATES, MB), # the rates of change
                    AUX) # additional output per time step
    return(outList)
  })
}

# Set parameter values that influence the initial value of W
Volume <- 0.0001  # 0.1*10^-3 m3
Cini   <- 10.0    # initial CO2 concentration in g CO2 m-3

# Solve model numerically with specified parameter values and variation in parameter K
states_ode45_K1 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 1, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

states_ode45_K10 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 10, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

states_ode45_K.1 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 0.1, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

# Plot dynamics of C as function of K:
plot(states_ode45_K10[, "time"], states_ode45_K10[, "C"],
  type = "l", xlim = c(0, 200), ylim = c(0, 11),
  xlab = "Time", ylab = "Conc. CO2", col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "C"], col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "C"], col = "red")

# add legend without a box
legend("topright", legend = c("K = 10", "K = 1", "K = 0.1"),
       col = c(" black", "blue", "red"), lty = c(1, 1, 1), bty = "n")

When $K$ is high (e.g. 10), photosynthesis is relatively slow which can be seen from the slow decrease in CO2 over time but also from the equation of $P$ ($P=...$)– larger $K$ leads to a small value of $C/(C+K)$, in turn leading to a smaller P. The smaller the K, the faster the photosynthesis which can be seen from the faster drop in CO2 over time but also from the equation of P – the closer $K$ gets to 0, the larger the value of $(C/(C+K)$, in turn leading to a larger $P$.

Negative concentrations are not physically possible. However, in the model it is possible! You may try this by setting the parameter $C_{ini}$ to a negative value.

When K equals 0, K drops out of the ratio ($C/(C+K)$) in the photosynthesis equation, and this ratio becomes 1 since C is in both the numerator and denominator, making P independent of the CO2 concentration in the chamber. One exception exist, which is when C is 0, then the ratio cannot be determined and will produce an error due to the division by zero. In other words, the photosynthesis rate P equals the maximum rate of photosynthesis $v_m$ when $K = 0$.

states_ode45_K0 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(1, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(Volume = Volume, Cini = Cini, K = 0), 
  method = "ode45"
)

plot(x = states_ode45_K1[, "time"], y = states_ode45_K1[, "C"],
  type = "l", xlim = c(0, 200), ylim = c(0, 11),
  xlab = "Time", ylab = "Conc. CO2", col = "blue")
lines(x = states_ode45_K10[, "time"], y =states_ode45_K10[, "C"], col = "black")
lines(x = states_ode45_K.1[, "time"], y =states_ode45_K.1[, "C"], col = "red")
lines(x = states_ode45_K0[, "time"], y =states_ode45_K0[, "C"], col = "green3")

# add legend without a box
legend("topright", legend = c("K = 10", "K = 1", "K = 0.1", "K = 0"),
  col = c("black", "blue", "red", "green3"), lty = 1, bty = "n")

The implementation of CUVET2:

# Define a function that returns a named vector with initial state values
InitCUVET2 <- function(W = 0){
  # Gather function arguments in a named vector
  y <- c(W = W)
  
  # Return the named vector
  return(y)
}

# Define a function that returns a named vector with parameter values
ParmsCUVET2 <- function(Volume = 0.0001,  # volume of the chamber, m3
                        Area   = 0.01,    # area of leaves, m2
                        Flow   = 0.00001, # through flow of air, m3/s
                        Cin    = 1.0,     # CO2 concentration of the air flowing in (g CO2 m-3)
                        Vm     = 0.001,   # Maximum rate of photosynthesis VM in g of CO2 m-2 s-1
                        K      = 1.0,     # K is the M-M constant of photosynthesis
                                          # with respect to CO2 in same units as C
                        Cini   = 10.00){  # initial CO2 concentration in g CO2 m-3
  # Gather function arguments in a named vector
  y <- c(Volume = Volume, 
         Area = Area, 
         Flow = Flow, 
         Cin = Cin, 
         Vm = Vm, 
         K = K, 
         Cini = Cini)
  
  # Return named vector with values
  return(y)
}

# Define the function returning the rates of change of the state variables with respect to time
RatesCUVET2 <- function(t, y, parms){
  with(as.list(c(y, parms)),{
    # Compute auxiliary/helper functions
    
    # Concentration of CO2 inside the chamber
    C <- W / Volume  
    # Net rate of photosynthesis (no respiration assumed)
    P <- Vm * C / (K + C)
    
    CIN  <- Cin * Flow # Inflow
    COUT <- C * Flow   # Outflow

    # Compute rate equations
    
    # Net rate of change of amount of CO2 in chamber
    dWdt <- CIN - COUT - P * Area        

    # Gather all rates of change in a vector
    RATES <- c(dWdt)
    
    # Optional: get in/out flow used to compute mass balances (or set MB <- NULL)
    MB <- NULL
    
    # Optional: gather auxiliary variables (as named vector) that should be returned
    AUX <- c(C = C, P = P, CIN = CIN, COUT = COUT)
    
    # Return result as a list
    outList <- list(c(RATES, MB), # first element: the rates of change
                    AUX) # additional output per time step
    return(outList)
  })
}

Solving the model for different values of $K$:

# Set parameter values that influence the initial value of W
Volume <- 0.0001 # m3
Cini   <- 0.0    # initial CO2 concentration in g CO2 m-3

# Solve model numerically with specified parameter values and variation in parameter K
states_ode45_K1 <- ode(
  y = InitCUVET2(W = Volume * Cini),
  times = seq(0, 200, by = 1),
  func = RatesCUVET2,
  parms = ParmsCUVET2(K = 1,
                      Volume = Volume,
                      Cini = Cini), # all others at their default values
  method = "ode45")

states_ode45_K10 <- ode(
  y = InitCUVET2(W = Volume * Cini),
  times = seq(0, 200, by = 1),
  func = RatesCUVET2,
  parms = ParmsCUVET2(K = 10,
                      Volume = Volume,
                      Cini = Cini), # all others at their default values),
  method = "ode45")

states_ode45_K.1 <- ode(
  y = InitCUVET2(W = Volume * Cini),
  times = seq(0, 200, by = 1),
  func = RatesCUVET2,
  parms = ParmsCUVET2(K = 0.1,
                      Volume = Volume,
                      Cini = Cini), # all others at their default values),
  method = "ode45")

# Plot dynamics of C as function of K:
plot(states_ode45_K10[, "time"], states_ode45_K10[, "C"],
  type = "l", xlim = c(0, 200), ylim = c(0, 1.2),
  xlab = "Time", ylab = "Conc. CO2", col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "C"], col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "C"], col = "red")

# add legend without a box
legend("topright", legend = c("K = 10", "K = 1", "K = 0.1"),
       col = c(" black", "blue", "red"), lty = c(1, 1, 1), bty = "n")

With a larger $K$ value, the photosynthetic rate becomes smaller, hence the uptake rate of CO₂ from the air by the leaf in the chamber is smaller. When the system is in a (dynamic) equilibrium the resulting concentrations of CO₂ will therefore be higher.

First, define a function (EquilCUVET2) that returns the equilibria (given parameter values as input):

EquilCUVET2 <- function(param){
  with(as.list(param),{
    # algebraic solution of the steady state equation (C-equilibrium)
    B   <- K - Cin + Vm * Area / Flow
    CEQUIL <- -B/2.0 + sqrt((B/2.0)^2 + Cin * K)
    return(CEQUIL)
  })
}

Then, compute the equilibria for the K values 0.1, 1 and 10:

Volume <- 0.0001 # m3
Cini   <- 0.0      # initial CO2 concentration in g CO2 m-3
EqK0.1 <- EquilCUVET2(ParmsCUVET2(K = 0.1, Volume = Volume, Cini = Cini))
EqK1.0 <- EquilCUVET2(ParmsCUVET2(K = 1,   Volume = Volume, Cini = Cini))
EqK10  <- EquilCUVET2(ParmsCUVET2(K = 10,  Volume = Volume, Cini = Cini))

Lastly, plot the results

par(mfrow = c(2, 1), # 2 rows, 1 column
    mar = c(2, 4, 0.1, 0.1)) # with margins = c(bottom, left, top, right)

plot(states_ode45_K10[, "time"], states_ode45_K10[, "C"],
  type = "l", xlim = c(0, 50), ylim = c(0, 1.2), 
  xlab = "Time", ylab = "Conc. CO2", lty = 1, lwd = 1, col = "black")
lines(c(-1E6, 1E6), c(EqK10, EqK10), lty = 2, lwd = 2, col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "C"], lty = 1, lwd = 1, col = "blue")
lines(c(-1E6, 1E6), c(EqK1.0, EqK1.0), lty = 2, lwd = 2, col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "C"], lty = 1, lwd = 1, col = "red")
lines(c(-1E6, 1E6), c(EqK0.1, EqK0.1), lty = 2, lwd = 2, col = "red")

legend("topleft", legend = c("K=10", "K=1.0", "K=0.1"), 
       col = c("black", "blue", "red"), lty = 1, bty = "n")

par(mar = c(4, 4, 0.1, 0.1))
plot(states_ode45_K10[, "time"], states_ode45_K10[, "P"],
  type = "l", xlim = c(0, 50), ylim = c(0, 8e-4), 
  xlab = "Time", ylab = "Rate of photosynthesis", col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "P"], col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "P"], col = "red")

legend("topleft", legend = c("K=10", "K=1.0", "K=0.1"), 
       col = c("black", "blue", "red"), lty = 1, bty = "n")

Eventually, the dynamic simulations result in the same equilibrium as the analytic result, providing confidence in the result. Note that the system reaches the equilibrium quicker when $K$ is 0.1 when compared to a $K$ of 10.

The time coefficient of depletion can be estimated in various ways, depending on the precision of the definition. When assuming that the leaf photosynthesizes is equal to the maximum rate Vm all the time, irrespective of the CO₂ concentration, we can estimate TC as follows: The initial amount of CO₂ equals $Volume * C_{ini}$ or 1 * 10^-3 g. The rate of depletion then equals $v_m * Area$ or 1 * 10^-5 g s^-1. The time to reach the equilibrium is then calculated from the ratio: TC = 100 s.

However, dW/dt is not constant and TC varies because $W$ decreases over time. To account for this, with $TC = \frac{W}{dW/dt}$, we rewrite dW/dt as a function of $W$:

$$P × Area$$

or

$$\frac{Area × v_m × C}{(C + K)}$$

or

$$\frac{Area × v_m × (W/Volume)}{(W/Volume ) + K}$$

if we multiply by 1, or more precisely by VOLUME/VOLUME, we get:

$$\frac{Area × v_m × W}{W + K × Volume}$$

since:

$$\frac{W}{dW/dt} = \frac{W}{(Area × v_m × W) / (W + K × Volume)}$$

and thus:

$$TC = \frac{W}{(Area × v_m × W) / (W + K × Volume)}$$

where W cancels:

$$TC = \frac{1}{(Area × v_m × 1) / (W + K × Volume)}$$

or,

$$TC = \frac{W + K × Volume}{Area × v_m}$$

Remarkably, the value of the time coefficient ($\tau$) decreases as W decreases, instead of being constant.

The upper asymptotic level for $\tau$, so where $K$ is very small, is $W/(Area * v_m)$ which is identical to our earlier calculation, because $W = Volume * C_{ini}$. The lower asymptotic limit for $\tau$, so where $W$ has become very small, is $(K * Volume) / (Area * v_m)$. For the smallest value of $K$ (=0.1) this is equal to 1 second. We may conclude that, although the model is quite simple, the determination of the time coefficient of the model and the appropriate time step is not easy and depends on parameter values. This is especially important to know when we would use Euler as integration method.

Added parts to function RatesCUVET1 as already defined above in exercise 5.1a):

# Compute the time coefficient
TC <- (W + K * Volume ) / (Area * Vm)

# Gather auxiliary information
AUX <- c(C = C, P = P, TC = TC)

Full function definition with these lines added:

# Define the function returning the rates of change of the state variables with respect to time
RatesCUVET1 <- function(t, y, parms){
  with(as.list(c(y, parms)),{
    # Compute auxiliary/helper functions

    # Concentration of CO2 inside the chamber
    C   <- W / Volume  
    # Net rate of photosynthesis (no respiration assumed)
    P <- Vm * C / (K + C)
    
    # Compute rate equations
    
    # Net rate of change of amount of CO2 in chamber
    dWdt <- - P * Area
    
    # Gather all rates of change in a vector
    # - the rates should be in the same order as the states (as specified in 'y')
    # - it can be a named vector, but does not need to be
    RATES <- c(dWdt = dWdt)
    
    # Optional: get in/out flow used to compute mass balances (or set MB <- NULL)
    # not included here (thus here use MB <- NULL), see template 3
    MB <- NULL

    # Compute the time coefficient
    TC <- (W + K * Volume ) / (Area * Vm)

    # Optional: gather auxiliary variables that should be returned (or set AUX <- NULL)
    # - this should be a named vector or list!
    AUX <- c(C = C, P = P, TC = TC)
    
    # Return result as a list
    # - the first element is a vector with the rates of change (in the same order as 'y')
    # - all other elements are (optional) extra output, which should be named
    outList <- list(c(RATES, MB), # the rates of change
                    AUX) # additional output per time step
    return(outList)
  })
}

We can plot the TC after solving the model in the same way we did in exercise 5.1b):

Volume <- 0.0001 # m3
Cini   <- 10.0   # initial CO2 concentration in g CO2 m-3

states_ode45_K1 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 1, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

states_ode45_K10 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 10, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

states_ode45_K.1 <- ode(
  y = InitCUVET1(W = Volume * Cini),
  times = seq(0, 2000, by = 1),
  func = RatesCUVET1,
  parms = ParmsCUVET1(K = 0.1, 
                      Volume = Volume, 
                      Cini = Cini), # all others at their default values
  method = "ode45")

Now we can plot it:

par(mfrow = c(2, 2), mar = c(4, 4, 1, 1)) # with mar=c(bottom, left, top, right)

plot(states_ode45_K10[, "time"], states_ode45_K10[, "C"],
  type = "l", xlim = c(0, 200), ylim = c(0, Cini), xlab = "Time", ylab = "CO2, g/m3", col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "C"], col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "C"], col = "red")

plot(states_ode45_K10[, "C"], states_ode45_K10[, "P"],
  type = "l", xlim = c(0, 10), ylim = c(0, 10E-4),
  xlab = "Conc. CO2", ylab = "Rate of photosynthesis", lty = 1, lwd = 1, col = "black")
lines(states_ode45_K1[, "C"], states_ode45_K1[, "P"], lty = 1, lwd = 1, col = "blue")
lines(states_ode45_K.1[, "C"], states_ode45_K.1[, "P"], lty = 1, lwd = 1, col = "red")

plot(states_ode45_K10[, "time"], states_ode45_K10[, "TC"],
  type = "l", xlim = c(0, 200), ylim = c(0, 200),
  xlab = "Time", ylab = "Time coefficient", lty = 1, lwd = 1, col = "black")
lines(states_ode45_K1[, "time"], states_ode45_K1[, "TC"], lty = 1, lwd = 1, col = "blue")
lines(states_ode45_K.1[, "time"], states_ode45_K.1[, "TC"], lty = 1, lwd = 1, col = "red")

# add legend without a box
legend("right", legend = c("K = 10", "K = 1", "K = 0.1"), 
       col = c(" black", "blue", "red"), lty = c(1, 1, 1), bty = "n")

plot(states_ode45_K10[, "W"], states_ode45_K10[, "TC"],
  type = "l", xlim = c(0, Volume * Cini), ylim = c(0, 200),
  xlab = "W", ylab = "Time coefficient", lty = 1, lwd = 1, col = "black")
lines(states_ode45_K1[, "W"], states_ode45_K1[, "TC"], lty = 1, lwd = 1, col = "blue")
lines(states_ode45_K.1[, "W"], states_ode45_K.1[, "TC"], lty = 1, lwd = 1, col = "red")

The plot of chamber concentration $C$ versus time, which shows a decreasing pattern with the abscissa (“x-axis”) as asymptote. The smaller the $K$-value, the steeper the curve. For a very small value of $K$, an almost linear decrease occurs until the CO₂ is practically depleted. The plot of $P$ versus $C$ shows the three hyperbolic Michaelis-Menten relationships, as they evolve over the time of the simulation period. There is a numerical problem connected with the fact that numerically negative values of $C$ may occur, even though this is physically not possible. In the negative range of $C$ the hyperbolic relationship has a vertical asymptote at $C=-K$. For even more negative values of $C$, the function value becomes positive again other “leg” of the hyperbola).