Land C cycle is a 1-box model. Steady-state in 1958. The sensitivity of the C-input flux to CO2 follows a logarithmic relationship to the CO2-change relative to 1958. A linear multiplication factor is applied to the sensitivity. Then, diagnose the sink flux in years 1980-2012.

1-box

## //////////////////////////////////////////////////
## 1-BOX MODEL
## --------------------------------------------------
onebox <- function( c_influx, tau, cpool0, return_cpool = TRUE) {
  ## ------------------------------------------------
  ## c_influx:  flux into pool (GtC a-1)
  ## cpool0:    initial pool size (GtC)
  ## tau:       turnover (residence) time (a)
  ## ------------------------------------------------
  tauisvec <- ifelse(length(tau)==length(c_influx), TRUE, FALSE)

  ## determine integration length (number of time steps) from length of 'c_influx'
  len <- length(c_influx)

  ## initialise output variable (time series of pool size)
  out_cpool <- rep( NA, len )
  out_outflux <- rep(NA, len)
  
  ## copy initial pool size to first element in output time series
  cpool <- cpool0

  ## integrate over each time step (this is an implementation of the differential equation)
  for (yr in seq(len) ) {

    ## copy current pool size to output time series
    out_cpool[yr] <- cpool
    
    ## outflux
    outflux <- 1/ifelse(tauisvec, tau[yr], tau) * cpool

    ## update pool size with input and decay
    cpool <- cpool + c_influx[yr] - outflux
    out_outflux[yr] <- outflux

  }

  ## function return value is a vector containing the output time series
  if (return_cpool){
    return( out_cpool )
  } else {
    return( out_outflux )
  }

}
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

# assume steady state in 1858 with C-input flux of 60, then response to CO2 change until 2012, then constant CO2
df_co2 <- read.table(paste0(here::here(), "/lab_bern/co2_monthly_maunaloa.txt"),
                      skip=70,
                      header=TRUE,
                      na.strings='-99.99'
                      ) |> 
  as_tibble() |> 
  dplyr::filter(mo == 1)  # take only january - 1 value per yer

get_sink <- function(sens, df_co2, years){
  
  c_influx_init <- 60
  d_influx <- sens * log(df_co2$co2_int/df_co2$co2_int[1])
  
  c_influx  <- c(rep( c_influx_init, 99 ),    # 100 years at constant C influx
                 c_influx_init + d_influx,    # transient response to CO2
                 rep(c_influx_init + d_influx[nrow(df_co2)], 300)  # 300 years stabilisation
                 )
  
  # run 1-box model
  df <- tibble(year = seq(length(c_influx)) + 1859 - 1, 
               c_influx = c_influx) %>%
    dplyr::mutate(c_pool = onebox(c_influx = c_influx, tau = 50, cpool0 = 3000)) |> 
    dplyr::mutate(dc = c_pool - lag(c_pool, default = 3000))
  
  # diagnose sink strength 1980-2012
  dc <- df |> 
    filter(year %in% years) |> 
    summarise(dc = mean(dc)) |> 
    pull(dc)
  
  return(dc)
}

vec_sens <- seq(from = 50, to = 150, by = 10)
vec_sink <- purrr::map_dbl(as.list(vec_sens),
           ~get_sink(., df_co2, years = 1980:2012))

tibble(sensitivity = vec_sens, sink = vec_sink) |> 
  ggplot(aes(sensitivity, sink)) +
  geom_point(color = "red", size = 3)