Data

library(dplyr)
library(tidyr)
library(ggplot2)
library(rsofun)
library(here)
library(lubridate)
library(FluxDataKit)

Read nice file from FluxDataKit. This contains meteorological data, remotely sensed LAI and fAPAR, and observations-based GPP from an extended set of FLUXNET sites.

path_forcingfile <- "~/data/FluxDataKit/v3/rsofun_driver_data_v3.rds"
driver <- readRDS(path_forcingfile)

Get data from all sites except:

sites <- FluxDataKit::fdk_site_info |>
  dplyr::filter(!(sitename %in% c("MX-Tes", "US-KS3"))) |>  # failed sites
  dplyr::filter(!(igbp_land_use %in% c("CRO", "WET"))) |> 
  left_join(
    fdk_site_fullyearsequence,
    by = "sitename"
  ) |> 
  filter(!drop) |> 
  filter(nyears >= 5)

Cut driver data to what’s filtered above.

driver_forcing_new <- driver |>
  select(sitename, forcing) |> 
  filter(sitename %in% sites$sitename) |> 
  left_join(
    sites |> 
      select(
        sitename, 
        year_start = year_start_fullyearsequence, 
        year_end = year_end_fullyearsequence),
    by = join_by(sitename)
  ) |> 
  unnest(forcing) |> 
  mutate(year = lubridate::year(date)) |> 
  mutate(keep = ifelse(year >= year_start & year <= year_end, TRUE, FALSE)) |> 
  dplyr::select(-year_start, -year_end, -year, -keep) |> 
  group_by(sitename) |> 
  nest()
  
driver <- driver |> 
  filter(sitename %in% sites$sitename)

driver$forcing <- driver_forcing_new$data

Model runs

Run the P-model (rsofun package) assuming 400 ppm for all sites and dates.

params_modl <- list(
  kphio              = 0.04998,
  kphio_par_a        = 0.0,
  kphio_par_b        = 1.0,
  soilm_thetastar    = 0.6 * 240,
  soilm_betao        = 0.0,
  beta_unitcostratio = 146.0,
  rd_to_vcmax        = 0.014,
  tau_acclim         = 30.0,
  kc_jmax            = 0.41
)

# run the model
output <- rsofun::runread_pmodel_f(
  driver |> 
    mutate(forcing = purrr::map(forcing, ~mutate(., co2 = 400))),
  par = params_modl
)

# some simulations failed due to missing values in the forcing. If that happens,
# the number of rows in the output data is 1.
output <- output |>
  mutate(len = purrr::map_int(data, ~nrow(.))) |>
  filter(len >= 1) |> # at least 3 years data
  select(-len)

Model run with varying CO2 as observed.

# run the model
output_varco2 <- rsofun::runread_pmodel_f(
  driver,
  par = params_modl
)
## Warning: There were 43 warnings in `dplyr::mutate()`.
## The first warning was:
## ℹ In argument: `data = purrr::pmap(., run_pmodel_f_bysite, params_modl = par,
##   makecheck = makecheck)`.
## Caused by warning in `FUN()`:
## ! Error: Missing value in temp for BE-Bra
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 42 remaining warnings.
# some simulations failed due to missing values in the forcing. If that happens,
# the number of rows in the output data is 1.
output_varco2 <- output_varco2 |>
  mutate(len = purrr::map_int(data, ~nrow(.))) |>
  filter(len >= 1) |> # at least 3 years data
  select(-len)

This yields data and model outputs for 216 sites. Let’s look for trends…