GPP trend and CO2?

Data

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

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/rsofun_driver_data_clean.rds"
driver <- readRDS(path_forcingfile)

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

Model run

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)

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

Trends in observed GPP

Annual mean

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

# fit linear regression to annual sums for each site and get slope
tmp <- driver |>
  mutate(adf = purrr::map(forcing, ~get_annual(.))) |> 
  mutate(linmod = purrr::map(adf, ~lm(gpp ~ year, data = .))) |> 
  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2]))

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

tmp |> 
  ggplot(aes(slope, after_stat(count))) +
  geom_histogram(fill = "grey", color = "black", bins = 50) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(tmp$slope), color = "red")

Main growing season

Mean across pper 33% quantile w.r.t. daily GPP values

get_annual <- function(df){
  threshold <- quantile(df$gpp, probs = 0.66, na.rm = TRUE)
  df |>
    filter(gpp > threshold) |> 
    mutate(year = lubridate::year(date)) |>
    group_by(year) |>
    summarise(gpp = mean(gpp))
}

# fit linear regression to annual sums for each site and get slope
tmp <- driver |>
  mutate(adf = purrr::map(forcing, ~get_annual(.))) |> 
  mutate(linmod = purrr::map(adf, ~lm(gpp ~ year, data = .))) |> 
  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2]))

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

tmp |> 
  ggplot(aes(slope, after_stat(count))) +
  geom_histogram(fill = "grey", color = "black", bins = 50) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(tmp$slope), color = "red")

De-seasonalised daily values

TBC

Trends in LUE

Annual mean

get_annual <- function(df){
  df |>
    mutate(year = lubridate::year(date),
           lue = gpp / (fapar * ppfd)) |>
    group_by(year) |>
    summarise(lue = mean(lue))
}

# fit linear regression to annual sums for each site and get slope
tmp <- driver |>
  filter(sitename != "SE-Svb") |> 
  mutate(adf = purrr::map(forcing, ~get_annual(.))) |> 
  mutate(linmod = purrr::map(adf, ~try(lm(lue ~ year, data = .)))) |> 
  mutate(error = purrr::map_chr(linmod, ~class(.))) |> 
  filter(error != "try-error") |> 
  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2])) |> 
  filter(slope > -12000) # outlier

## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases
## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases
## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases
## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases
## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases
## Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
##   0 (non-NA) cases

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

tmp |> 
  ggplot(aes(slope, after_stat(count))) +
  geom_histogram(fill = "grey", color = "black", bins = 50) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(tmp$slope), color = "red")

Main growing season

TBC

De-seasonalised daily values

TBC

Trends in P-model residuals

Combine outputs and forcing.

tmp <- output |> 
  select(sitename, data) |> 
  unnest(data) |> 
  rename(gpp_pmodel = gpp) |> 
  left_join(
    driver |> 
      select(sitename, forcing) |> 
      unnest(forcing),
    by = c("sitename", "date")
  ) |> 
  mutate(res = gpp_pmodel - gpp) |>
  group_by(sitename) |> 
  nest()

No need to do annual totals. Fit linear regression directly on residuals.

# fit linear regression to annual sums for each site and get slope
tmp2 <- tmp |>
  mutate(len = purrr::map_int(data, ~nrow(.))) |>
  filter(len >= 1000) |> # at least 3 years data
  select(-len) |> 
  mutate(linmod = purrr::map(data, ~lm(res ~ date, data = .))) |> 
  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2]))

Plot distribution of linear trends in residuals - no dominance of negative trends (expected since the model by design doesn’t account for the positive CO₂ effect on GPP).

tmp2 |> 
  ggplot(aes(slope, after_stat(count))) +
  geom_histogram(fill = "grey", color = "black", bins = 50) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(tmp2$slope), color = "red")

Main growing season

filter_gs <- function(df){
  threshold <- quantile(df$gpp, probs = 0.66, na.rm = TRUE)
  df |>
    mutate(res = ifelse(gpp > threshold, res, NA))
}

# fit linear regression to annual sums for each site and get slope
tmp2 <- tmp |>
  mutate(len = purrr::map_int(data, ~nrow(.))) |>
  filter(len >= 1000) |> # at least 3 years data
  select(-len) |> 
  mutate(data = purrr::map(data, ~filter_gs(.))) |> 
  mutate(linmod = purrr::map(data, ~lm(res ~ date, data = .))) |> 
  mutate(slope = purrr::map_dbl(linmod, ~coef(.)[2]))

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

tmp2 |> 
  ggplot(aes(slope, after_stat(count))) +
  geom_histogram(fill = "grey", color = "black", bins = 50) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_vline(xintercept = median(tmp2$slope), color = "red")

De-seasonalised daily values

TBC