GPP trend and CO2?

Data

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

# if(!require(devtools)){install.packages(devtools)}
# devtools::install_github("geco-bern/rsofun")
library(rsofun)

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_2/FluxDataKit/v3.4/zenodo_upload/rsofun_driver_data_v3.4.2.rds"
driver <- readRDS(path_forcingfile)

Get data from all sites except:

Croplands
Wetlands
Less than 5 years of sequential good-quality data

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_gpp) |> 
  filter(nyears_gpp >= 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 = start_gpp, 
        year_end = end_gpp),
    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.

# values from sofunCalVal reference run
params_modl <- list(
  kphio              = 0.046857840,
  kphio_par_a        = -0.001836235,
  kphio_par_b        = 19.866210640,
  soilm_thetastar    = 28.055517785,
  soilm_betao        = 0.006899699,
  beta_unitcostratio = 146.0,
  rd_to_vcmax        = 0.014,
  tau_acclim         = 30.0,
  kc_jmax            = 0.41
)

# run the model
output_constco2 <- 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_constco2 <- output_constco2 |>
  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
)

# 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…

Trends in GPP

get_annual <- function(df){
  df |>
    mutate(year = lubridate::year(date)) |>
    group_by(year) |>
    summarise(gpp = sum(gpp))
}

get_mean <- function(df){
  mean(df$gpp)
}

# fit linear regression to annual sums for each site and get slope
df_obs_ann <- driver |>
  mutate(ndata_gpp = purrr::map_int(forcing, ~sum(!is.na(.$gpp)))) |> 
  filter(ndata_gpp > 0) |> 
  mutate(adf = purrr::map(forcing, ~get_annual(.))) |> 
  mutate(linmod = purrr::map(adf, ~lm(gpp ~ year, data = .))) |> 
  mutate(slope_obs = purrr::map_dbl(linmod, ~coef(.)[2])) |> 
  mutate(mean_obs = purrr::map_dbl(adf, ~get_mean(.))) |> 
  mutate(slope_rel_obs = slope_obs / mean_obs) |> 
  
  # remove outlier
  filter(slope_obs > -150)

df_mod_ann <- output_varco2 |> 
  mutate(ndata_gpp = purrr::map_int(data, ~sum(!is.na(.$gpp)))) |> 
  filter(ndata_gpp > 0) |> 
  mutate(adf = purrr::map(data, ~get_annual(.))) |> 
  mutate(linmod = purrr::map(adf, ~lm(gpp ~ year, data = .))) |> 
  mutate(slope_mod = purrr::map_dbl(linmod, ~coef(.)[2])) |> 
  mutate(mean_mod = purrr::map_dbl(adf, ~get_mean(.))) |> 
  mutate(slope_rel_mod = slope_mod / mean_mod)

df_ann <- df_obs_ann |> 
  select(sitename, slope_obs, mean_obs, slope_rel_obs) |> 
  left_join(
    df_mod_ann |> 
      select(sitename, slope_mod, mean_mod, slope_rel_mod),
    by = "sitename"
  )

Absolute trends

Plot distribution of linear trends - no dominance of positive trends.

df_ann |> 
  ggplot() +
  geom_histogram(
    aes(slope_obs, after_stat(count), fill = "Observations"), 
    color = "black",
    alpha = 0.5,
    bins = 30
    ) +
  geom_histogram(
    aes(slope_mod, after_stat(count), fill = "P-model"), 
    color = "black", 
    alpha = 0.5,
    bins = 30
    ) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(df_ann$slope_obs, na.rm = TRUE), color = "black") +
  geom_vline(xintercept = median(df_ann$slope_mod, na.rm = TRUE), color = "red") +
  theme_classic() +
  labs(x = expression(paste("Absolute GPP trend (gC m"^-2, " yr"^-2, ")"))) +
  scale_fill_manual(
    name = "",
    values = c(
      "Observations" = "black",
      "P-model" = "red"
    )
  )

df_ann |> 
  ggplot(aes(slope_mod, slope_obs)) +
  geom_point() +
  labs(
    x = expression(paste("Modeled absolute GPP trend (gC m"^-2, " yr"^-2, ")")),
    y = expression(paste("Observed absolute GPP trend (gC m"^-2, " yr"^-2, ")"))
    ) +
  geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
  coord_equal() +
  theme_classic()

Relative trends