Climate

We retrieved all data from TerraClimate detailed below (Abatzoglou et al. 2018). This could be compared to site data, but with issues regarding site data standardization, representativeness, and validation.

ET: Actual evapotranspiration, mm
CWD: Climate water deficit, mm
Tmax: Maximum temperature, °C
Tmin: Minimum temperature, °C
PDSI: Palmer Drought Severity Index
PET: Potential evapotranpiration, mm
Pr: precipitation, mm
VPD: Vapour Pressure Deficit, kPa

Code

list.files("data/climate/", full.names = TRUE) %>%
  read_tsv() %>%
  group_by(site) %>%
  mutate(
    aet = aet * .1, def = def * .1, pdsi = pdsi * .01, pet = pet * .1,
    tmmn = tmmn * .1, tmmx = tmmx * .1, vpd = vpd * .01
  ) %>%
  write_tsv("outputs/climate.tsv")

For a quick exploration of each site climatology we first had a look to umbrothermal diagrams for the mean climate of each site, opposing for instance highly seasonal monsoon regimes from India un Uppangala to less seasonal regimes from Missiones in Argentina.

Code

read_tsv("outputs/climate.tsv") %>%
  select(-lon, -lat) %>%
  group_by(site, month) %>%
  summarise_all(mean) %>%
  ggplot(aes(x = month)) +
  geom_col(aes(y = pr / 30), fill = "grey", col = NA) +
  geom_line(aes(y = tmmn, col = "maximum")) +
  geom_line(aes(y = tmmx, col = "minimum")) +
  theme_bw() +
  xlab("") +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
  ) +
  scale_y_continuous(
    name = "Temperature (°C)",
    sec.axis = sec_axis(trans = ~ . * 30, name = "Precipitation (mm)")
  ) +
  theme(
    axis.title.y = element_text(color = "black"),
    axis.title.y.right = element_text(color = "grey")
  ) +
  scale_color_discrete("") +
  theme(legend.position = "bottom") +
  facet_wrap(~site)

Umbrothermal diagrams for the mean climate of each site.

Looking more in depth to each variable per site, plot and month better showed the climatology for each variable with relative low variation within site but large variation among sites.

Code

read_tsv("outputs/climate.tsv") %>%
  gather(variable, value, -site, -month, -plot, -lat, -lon) %>%
  mutate(var_long = recode(variable,
    "aet" = "Actual evapotranspiration [ mm ]",
    "def" = "Climate water deficit [ mm ]",
    "pdsi" = "Palmer Drought Severity Index",
    "pet" = "Potential evapotranspiration [ mm ]",
    "pr" = "Precipitation [ mm ]",
    "tmmn" = "Minimum temperature [ °C ]",
    "tmmx" = "Maximum temperature [ °C ]",
    "vpd" = "Vapor pressure deficit [ kPA ]"
  )) %>%
  ggplot(aes(month, value, col = site, group = paste(site, plot))) +
  geom_line() +
  facet_wrap(~var_long, scales = "free_y") +
  xlab("") +
  ylab("") +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
  ) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_color_discrete("")

Climate variables per site, plot and month.

To ease plot and sites comparisons, we derived yearly metrics of sum (evapotranspiration, climatic water deficit, precipitation), minimum (minimum temperature), or maximum (palmer drought index, maximum temperature, vapour pressure deficit) depending on the variables. This better revealed among sites differences with for instance highest CWD at Uppangala despite highest total precipitation.

Other indices could be derived such as widely used bioclimatic variables.

Code

read_tsv("outputs/climate.tsv") %>%
  group_by(site, plot) %>%
  summarise(
    aet = sum(aet), def = sum(def), pdsi = max(pdsi), pet = sum(pet),
    pr = sum(pr), tmax = max(tmmx), tmin = min(tmmn), vpd_max = max(vpd)
  ) %>%
  gather(variable, value, -site, -plot) %>%
  mutate(var_long = recode(variable,
    "aet" = "Evapotranspiration [ mm ]",
    "def" = "Climate water deficit [ mm ]",
    "pdsi" = "Palmer Drought Severity Index",
    "pdsi_max" = "Max PDSI",
    "pdsi_mean" = "Mean PDSI",
    "pet" = "Pot evapotranspiration [ mm ]",
    "pr" = "Precipitation [ mm ]",
    "tmin" = "Min temperature [ °C ]",
    "tmax" = "Max temperature [ °C ]",
    "vpd_max" = "Max VPD [ kPA ]"
  )) %>%
  ggplot(aes(site, value, fill = site, group = paste(site, plot))) +
  geom_col(position = "dodge") +
  facet_wrap(~var_long, scales = "free_x") +
  theme_bw() +
  coord_flip() +
  theme(legend.position = "none", axis.title = element_blank(),
        axis.text.y = element_text(size = 6))

The first ordination axis retained almost half of the variations of the 8 climate variables.

This is a very quick interpretation of the figure and should be explored more in depth.

Code

library(ggfortify)
data <- read_tsv("outputs/climate.tsv") %>%
  group_by(site, plot) %>%
  summarise(
    et = sum(aet), cwd = sum(def), pdsi = max(pdsi), pet = sum(pet),
    pr = sum(pr), tmax = max(tmmx), tmin = min(tmmn), vpd = max(vpd)
  ) %>%
  ungroup()
pca <- prcomp(data %>% select(-site, -plot), scale. = TRUE)
autoplot(pca,
  loadings = TRUE, loadings.label = TRUE,
  loadings.label.repel = TRUE,
  data = data, colour = "site"
) +
  theme_bw() +
  coord_equal() +
  scale_color_discrete("")

Yearly mean climate variables PCA per site and plot.