Setup

All parameters, site info, and forcing data is available for an example in the driver object (part of rsofun package)

  lm3ppa_p_model_drivers$forcing[[1]]
## # A tibble: 365 × 18
##    date        year   doy  hour    temp temp_s…¹    prec snow  vpd      rh  ppfd
##    <date>     <dbl> <dbl> <dbl>   <dbl>    <dbl>   <dbl> <lgl> <lgl> <dbl> <dbl>
##  1 2009-01-01  2009     1  11.5  0.384      3.39 1.66e-5 NA    NA     93.7  28.0
##  2 2009-01-02  2009     2  11.5 -1.64       2.55 2.32e-5 NA    NA     92.5  24.4
##  3 2009-01-03  2009     3  11.5 -2.51       1.95 3.71e-6 NA    NA     85.1  46.5
##  4 2009-01-04  2009     4  11.5 -1.82       2.20 1.30e-5 NA    NA     83.5  34.6
##  5 2009-01-05  2009     5  11.5 -1.34       2.23 2.23e-5 NA    NA     87.8  36.2
##  6 2009-01-06  2009     6  11.5 -0.450      2.75 2.19e-5 NA    NA     90.9  25.9
##  7 2009-01-07  2009     7  11.5  0.266      3.50 1.36e-5 NA    NA     89.7  34.8
##  8 2009-01-08  2009     8  11.5  0.504      3.52 1.13e-5 NA    NA     86.1  44.4
##  9 2009-01-09  2009     9  11.5  0.0869     3.49 1.86e-5 NA    NA     90.3  25.3
## 10 2009-01-10  2009    10  11.5 -0.404      3.44 1.25e-5 NA    NA     90.1  40.5
## # … with 355 more rows, 7 more variables: par <dbl>, patm <dbl>, wind <dbl>,
## #   ccov_int <lgl>, ccov <lgl>, co2 <dbl>, swc <dbl>, and abbreviated variable
## #   name ¹​temp_soil
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names

This contains the forcing time series data frame where the disturbance is to be defined as the fraction of aboveground biomass harvested (harv). Additional specifications of the disturbance forcing, which are more specific to the simulations done here, are hard-coded (see below).

Specify forcing

Disturbance regime

Model is first run to steady state. Then for another 100 years undisturbed (continued steady-state) until first disturbance. After first disturbance left undisturbed for 900 years to allow for recovery (in some simulations recovery may take long). Then a regime of repeated disturbance after simulation year 1000 to investigate disturbance-recovery (non-steady state) dynamics, first at low frequency for 1000 years (disturbance every 250 years), then at high frequency for 1000 years (disturbance every 25 years).

Disturbance is implemented by: - year 100 first disturbance applied, then recovery until year 1000 - year 1000 second disturbance, then every 250 years disturbed for 1000 years - … then every 25 years disturbed for another 1000 years

To be handled by model by (new) forcing time series as input (harv)

fharv <- 0.9
harv_vec <- rep(0, 999)
harv_vec[100] <- fharv
harv_vec <- c(harv_vec, rep(c(fharv, rep(0, 249)), 4), rep(c(fharv, rep(0, 24)), 40), 0)

df_harv <- tibble(year = seq(length(harv_vec)), harv = harv_vec)

df_harv %>% 
  ggplot(aes(year, harv)) +
  geom_line() +
  ylim(0, 1)

Experimental design

df_sims <- read_csv("simulations_disturb.csv")
## Rows: 13 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): Simulation name, C-N coupling, Disturbance, ABG removed from system...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
df_sims %>% knitr::kable()
Simulation name C-N coupling Disturbance ABG removed from system SOM removed Grass regrowing Soil temperature change N deposition level
sc0 No No 0% 0% No none control
sc1 No Yes 0% 0% No none control
sc2 No Yes 90% 0% No none control
sc3 No Yes 90% 0% Yes none control
sc4 No Yes 90% 0% No 5 deg C control
cn0 Yes No 0% 0% No none control
cn1 Yes Yes 0% 0% No none control
cn2 Yes Yes 90% 0% No none control
cn3 Yes Yes 90% 0% Yes none control
cn4 Yes Yes 90% 0% No 5 deg C control
cn5 Yes Yes 90% 50% No none control
cn6 Yes Yes 90% 0% No none Low
cn7 Yes Yes 90% 0% No none High

Create forcing objects

## get mean seasonal cycle and repeat this every year of all simulations
df_forcing <- lm3ppa_p_model_drivers$forcing[[1]] %>% 
  mutate(doy = lubridate::yday(date)) %>% 
  group_by(doy) %>% 
  summarise(across(is.numeric, mean))
## Warning: Predicate functions must be wrapped in `where()`.
## 
##   # Bad
##   data %>% select(is.numeric)
## 
##   # Good
##   data %>% select(where(is.numeric))
## 
## ℹ Please update your code.
## This message is displayed once per session.

Repeat mean seasonal cycle nyears times where nyears corresponds to the length of the harvest time series (rows in df_harv). The column year now signifies simulation year and goes from 1 to nyears.

Add harvest forcing to drivers.

nyears <- nrow(df_harv)

df_forcing <- df_forcing %>% 
  slice(rep(1:n(), nyears)) %>% 
  mutate(year = rep(1:nyears, each = 365))

Add harvest to forcing, assuming harvest on Jan 1st.

df_forcing_disturb <- df_forcing %>% 
  left_join(
    df_harv %>% 
      mutate(doy = 1),
    by = c("doy", "year")
  ) %>% 
  mutate(harv = ifelse(is.na(harv), 0, harv))
  
  # ## add pseudo-date, starting in year 2000
  # mutate(date = lubridate::ymd("0000-01-01") + lubridate::years(year-1) + lubridate::days(doy-1))

## for control simulation
df_forcing <- df_forcing %>%
  mutate(harv = 0)

Add N deposition as NOx and NHy.

df_forcing <- df_forcing %>% 
  mutate(nox = 0, nhy = 0)

df_forcing_disturb <- df_forcing_disturb %>% 
  mutate(nox = 0, nhy = 0)

Create new driver objects.

## control simulations without disturbance
lm3ppa_p_model_drivers_xx0 <- lm3ppa_p_model_drivers
lm3ppa_p_model_drivers_xx0$forcing[[1]] <- df_forcing
lm3ppa_p_model_drivers_xx0$params_siml[[1]]$firstyeartrend <- 0
lm3ppa_p_model_drivers_xx0$params_siml[[1]]$nyeartrend <- 3000

## simulations with disturbance
lm3ppa_p_model_drivers_xx1 <- lm3ppa_p_model_drivers
lm3ppa_p_model_drivers_xx1$forcing[[1]] <- df_forcing_disturb
lm3ppa_p_model_drivers_xx1$params_siml[[1]]$firstyeartrend <- 0
lm3ppa_p_model_drivers_xx1$params_siml[[1]]$nyeartrend <- 3000

## xxx debug
lm3ppa_p_model_drivers_xx1$forcing[[1]]$harv[1] <- 0

Model runs

sc1

filn <- "data/out_sc1.rds"
if (!file.exists(filn)){
  out_sc1 <- runread_lm3ppa_f(
       lm3ppa_p_model_drivers_xx1,
       makecheck = TRUE,
       parallel = FALSE
       )
  saveRDS(out_sc1, file = filn)
} else 
  out_sc1 <- readRDS(filn)

sc2

Export of dead biomass from system (not added to soil) is implemented by simply not calling the plant2soil() in the subroutine disturb() (file vegetation_lm3ppa.mod.f90). Comment it out and re-compile before running.

filn <- "data/out_sc2.rds"
if (!file.exists(filn)){
  out_sc2 <- runread_lm3ppa_f(
       lm3ppa_p_model_drivers_xx1,
       makecheck = TRUE,
       parallel = FALSE
       )
  saveRDS(out_sc2, filn)
} else 
  out_sc2 <- readRDS(filn)

sc3

Not done (for now).

sc4

The soil is heated compared to air temperatures if shading effect of trees and their canopy is removed. We can parametrise this as a function of the radiation reaching the ground, which is affected by top-of-canopy radiation and LAI. The function chosen here is to model \(\Delta T\) as a linear function of the solar radiation reaching the ground, so that it attains 5 deg C at 250 W m-2 (mean daily) and 0 at 0 W m-2.

The function chosen here is a sigmoid function with parameters set manually to get it into the right range (visually, see below).

lm3ppa_p_model_drivers_xx1$forcing[[1]] %>% 
  slice(1:365) %>% 
  ggplot(aes(doy, ppfd)) +
  geom_line()

lm3ppa_p_model_drivers_xx1$forcing[[1]] %>% 
  slice(1:365) %>% 
  ggplot(aes(doy, 5 * ppfd/250)) +
  geom_line()

filn <- "data/out_sc4.rds"
if (!file.exists(filn)){
  out_sc4 <- runread_lm3ppa_f(
       lm3ppa_p_model_drivers_xx1,
       makecheck = TRUE,
       parallel = FALSE
       )
  saveRDS(out_sc4, filn)
} else {
  out_sc4 <- readRDS(filn)
}

Analysis

Take annual tile-level outputs (for shorthand).

out_sc1_ann <- out_sc1$data[[1]]$output_annual_tile
out_sc2_ann <- out_sc2$data[[1]]$output_annual_tile
out_sc4_ann <- out_sc4$data[[1]]$output_annual_tile

Plant C

# model output includes the spinup. Remove it for plotting and overwrite years.
out_sc1_ann <- out_sc1_ann %>%
  slice((lm3ppa_p_model_drivers_xx1$params_siml[[1]]$spinupyears + 1):nrow(out_sc1_ann)) %>% 
  mutate(year = 1:nyears)
out_sc2_ann <- out_sc2_ann %>%
  slice((lm3ppa_p_model_drivers_xx1$params_siml[[1]]$spinupyears + 1):nrow(out_sc2_ann)) %>% 
  mutate(year = 1:nyears)
out_sc4_ann <- out_sc4_ann %>%
  slice((lm3ppa_p_model_drivers_xx1$params_siml[[1]]$spinupyears + 1):nrow(out_sc4_ann)) %>% 
  mutate(year = 1:nyears)

gg <- out_sc1_ann %>% 
  mutate(simulation = "sc1") %>% 
  bind_rows(
    out_sc2_ann %>% 
      mutate(simulation = "sc2")
  ) %>% 
  bind_rows(
    out_sc4_ann %>% 
      mutate(simulation = "sc4")
  ) %>% 
  ggplot() +
  geom_line(aes(x = year, y = plantC, color = simulation)) +
  theme_classic() +
  geom_vline(xintercept = df_harv %>% filter(harv > 0) %>% pull(year), alpha = 0.3) +
  labs(x = "Year", y = "Plant C") +
  ylim(0, 15)

gg

Soil C

gg <- out_sc1_ann %>% 
  mutate(simulation = "sc1") %>% 
  bind_rows(
    out_sc2_ann %>% 
      mutate(simulation = "sc2")
  ) %>% 
  bind_rows(
    out_sc4_ann %>% 
      mutate(simulation = "sc4")
  ) %>% 
  ggplot() +
  geom_line(aes(x = year, y = fastSOM + SlowSOM, color = simulation)) +
  theme_classic() +
  geom_vline(xintercept = df_harv %>% filter(harv > 0) %>% pull(year), alpha = 0.3) +
  labs(x = "Year", y = "Soil C") +
  ylim(0, 100)

gg