Code for this is available on my github.

Approach

Gao et al. (2014) suggest that the plant rooting is adapted to the cumulative water deficit during the dry spells They adopt the Mass Curve Technique to derive it. A Gumbel distribution is fit to the accumulated water deficits during dry spells and allows for an estimate of the deficit with a given return period. The method requires precipitation (\(P\)), runoff (\(Q\)), potential evapotranspiration (\(E_p\)) and green vegetation cover (\(f_v\)) to be specified for each time step (here daily) over multiple years.

The approach implemented here (function mct()) considers temporal variations in the demand, while Gao et al. calculated a mean annual and mean dry season demand. Here, PET is calculated based on the Priestly-Taylor equation as opposed to the Hargreaves equation used by Gao et al. Limitations of the method described here include that lateral runon and runoff and delayed water inflow by snow melt are ignored, and PET is assumed to drive the demand but is not affected by other plant adaptations to dry conditions (reduction of stomatal conductance, internal water storage). Effects of phenological changes on water demand during droughts are accounted for by \(f_v\).

The steps for the method are:

  1. Identify events where the water deficit (\(f_v E_p - P\)) is accumulating.
  2. Fit a gumbel distribution to the largest \(N\) events.
  3. Extract the estimated water deficit for a given return period \(T\).

Identify events

Let’s try this out for the FLUXNET site ‘FR-Pue’, for which whe have these variables.

# ## get data
# load("~/eval_pmodel/data/v3/out_eval_FULL.Rdata")
# df <- obs_eval_NT$ddf %>% 
#   as_tibble() %>% 
#   filter(sitename=="FR-Pue") %>% 
#   select(date, fv=fapar, prec=prec_fluxnet2015)
# 
# ## add PET from the model output
# load("~/eval_pmodel/data/mod_FULL.Rdata")
# df <- df %>% 
#   left_join( select( mod_FULL$daily$`FR-Pue`, date, pet ), by="date" ) %>% 
#   cutna_headtail_df("fv") %>% 
#   mutate(fv = rbeni::myapprox(fv), pet = rbeni::myapprox(pet), prec=ifelse(is.na(prec),0,prec))
# save(df, file = "data/df.Rdata")

load("data/df.Rdata")
df %>% 
  ggplot(aes(date, fv*pet)) +
  geom_line()

Plot the demand curve, given by cumulative \(f_v E_p\).

df <- df %>% 
  mutate(demand_cum = cumsum(fv * pet),
         supply_cum = cumsum(prec)) %>% 
  mutate(bal = prec - fv*pet) %>% 
  mutate(bal = myapprox(bal)) %>% 
  mutate(bal_cum = cumsum(bal))

# df %>% 
#   ggplot(aes(date, bal)) +
#   geom_line()

Get events of consecutive deficit using the mct() function.

source("R/mct.R")
out1 <- mct(df)
out2 <- mct(df, method = "threshbal", thresh_deficit=0.5)

Plot the cumulative deficit and rain events.

ggplot() +
  geom_rect(
    data=out2$inst, 
    aes(xmin=date_start, xmax=date_end, ymin=-99, ymax=99999), 
    fill=rgb(0,0,0,0.3), 
    color=NA) +
  # geom_line(data = out1$df, aes(date, -bal_cum), size = 0.3, color="royalblue") +
  geom_line(data = out1$df, aes(date, prec), size = 0.3, color="royalblue") +
  geom_line(data = out1$df, aes(date, deficit)) +
  geom_line(data = out2$df, aes(date, deficit), color="tomato") +
  coord_cartesian(ylim=c(0, 450), xlim = c(ymd("2002-01-01"), ymd("2010-12-01"))) +
  theme_classic()

As an additional illustration, plot cumulative variables.

ggplot() +
  geom_rect(
    data=out2$inst, 
    aes(xmin=date_start, xmax=date_end, ymin=-99, ymax=99999), 
    fill=rgb(0,0,0,0.3), 
    color=NA) + 
  geom_line(data=df, aes(date, demand_cum)) +
  geom_line(data=df, aes(date, supply_cum, color="supply_cum"), color="tomato") +
  labs(y=expression(integral(f[v] ~ E[p])), x="Date") +
  coord_cartesian(ylim=c(0, 10000)) +
  theme_classic()

Plot the distribution of cumulative deficits

ggplot(out2$inst, aes(deficit)) +
  geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Fit a distribution

To estimate the probability of extreme values, we fit a Gumbel distribution following this link. A paper on the R package extRemes.

## Take only the N largest instances (deficits), where N is given by the number of years available in the data
nyears <- year(range(out2$df$date)[2]) - year(range(out2$df$date)[1]) + 1
vals <- out2$inst %>% 
  arrange(desc(deficit)) %>% 
  dplyr::slice(1:nyears) %>% 
  select(deficit) %>% 
  unlist() %>% 
  unname()

gumbi <- extRemes::fevd(x=vals, type="Gumbel", method="MLE")
summary(gumbi)
## 
## extRemes::fevd(x = vals, type = "Gumbel", method = "MLE")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  67.08301 
## 
## 
##  Estimated parameters:
##  location     scale 
## 204.46074  35.40913 
## 
##  Standard Error Estimates:
##  location     scale 
## 10.336948  7.946064 
## 
##  Estimated parameter covariance matrix.
##          location    scale
## location 106.8525 25.56470
## scale     25.5647 63.13993
## 
##  AIC = 138.166 
## 
##  BIC = 139.2959
# extract MLEs (these are needed for the remaining part of the analysis)
muG <- gumbi$results$par[1]
sigmaG <- gumbi$results$par[2]

## don't know which package these are coming from (e1071?)
probplot(values=vals, model=gumbi, varname="Deficit (mm)", alpha=1-0.95, dist="gumbel")

#cprobplot(values=vals, model=gumbi, varname="Deficit (mm)", alpha=1-0.95, dist="gumbel")

QQplot(values=vals, mu=muG, sigma=sigmaG, dist="gumbel")

PPplot(values=vals, mu=muG, sigma=sigmaG, dist="gumbel")

# all plots ("primary")
# extRemes::plot.fevd(gumbi)

# only return period plot
plot.fevd.mle(gumbi, 
      # type = c("primary", "probprob", "qq", "qq2", "Zplot", "hist", "density", "rl", "trace"),
      type = c("rl"),
      rperiods = c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800),
      a = 0, hist.args = NULL, density.args = NULL, d = NULL )

source("R/get_return_period.R")
df_test <- get_return_period(gumbi)
## Warning in log(1 - 1/return_periods): NaNs produced
with(df_test, plot(trans_period, return_values, pch=16, col="red"))

Extract value for return periods

Get tranformed variate of return period \(T\) as described in Gao et al. as \[ y = - \ln ( -\ln (1-1/T) ) \]

## get return levels for a given vector of return periods
return_period <- c(2, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 200, 250, 300, 500, 800)
return_level <- extRemes::return.level(
  gumbi, 
  return.period = return_period
  )
df_return <- tibble( 
  return_period = return_period, 
  return_level = unname(c(return_level)), 
  trans_period = -log( -log(1 - 1/return_period)) )

df_return %>% 
  ggplot(aes(trans_period, return_level)) +
  geom_point()

Nice function

source("R/get_data_fluxnet2015_mct.R")
source("R/get_plantwhc_gao.R")
ddf <- get_data_fluxnet2015_mct("FR-Pue")
result <- get_plantwhc_mct(ddf)