Cumulative water deficit example

library(dplyr)
library(lubridate)
library(patchwork)
library(extRemes)
library(ggplot2)
library(cwd)

This demonstrates the workflow for determining cumulative water deficit (CWD) time series and fitting an extreme value distribution to annual maxima of the CWD time series.

Prepare data

Read data from the file contained in this repository.

df <- cwd:::df_CH_LAE |>
  # convert to from kPA to Pa (SI units are used for inputs to functions in the package)
  mutate(PA_F = 1e3 * PA_F)

Convert ET to mass units

Convert latent heat flux (W/m2) to evapotranspiration in mass units (mm/d).

# tested: identical results are obtained with:
# bigleaf::LE.to.ET(LE_F_MDS, TA_F_MDS)* 60 * 60 * 24
le_to_et <- function(le, tc, patm){
  1000 * 60 * 60 * 24 * le / (cwd::calc_enthalpy_vap(tc) * cwd::calc_density_h2o(tc, patm))
}

df <- df |>
  mutate(et = le_to_et(LE_F_MDS, TA_F_MDS, PA_F))

Check annual totals.

adf <- df |>
  mutate(year = year(TIMESTAMP)) |>
  group_by(year) |>
  summarise(et = sum(et), prec = sum(P_F))

adf |>
  tidyr::pivot_longer(cols = c(et, prec), names_to = "Flux") |> 
  ggplot(aes(x = year, y = value, color = Flux)) +
  geom_line() +
  labs(y = "Flux (mm/yr)")

Each year, the annual precipitation is greater than ET. Hence, the water deficit will not continue accumulating over multiple years.

Simulate snow

Simulate snow accumulation and melt based on temperature and precipitation.

df <- df |>
  mutate(prec = ifelse(TA_F_MDS < 0, 0, P_F),
         snow = ifelse(TA_F_MDS < 0, P_F, 0)) |>
  cwd::simulate_snow(varnam_prec = "prec", varnam_snow = "snow", varnam_temp = "TA_F_MDS")

Visualise snow mass equivalent time series.

df |> 
  ggplot(aes(TIMESTAMP, snow_pool)) + 
  geom_line() +
  labs(x = "Date", y = "Snow mass equivalent (mm)")

This looks like it’s a lot of snow, actually. Maybe the melting rate is too slow.

Define the daily water balance as liquid water infiltrating into soil (taken as rain plus snow melt) minus evapotranspiration - both in mass units, or equivalently in mm/d.

df <- df |>
  mutate(wbal = liquid_to_soil - et)

Visualise it.

gg5 <- df |>
  ggplot(aes(TIMESTAMP, et)) +
  geom_line()

gg6 <- df |>
  ggplot(aes(TIMESTAMP, liquid_to_soil)) +
  geom_line()

gg7 <- df |>
  ggplot(aes(TIMESTAMP, wbal)) +
  geom_line()

gg5 / gg6 / gg7

Cumulative water deficit algorithm

Get CWD and events.

out_cwd <- cwd(
  df,
  varname_wbal = "wbal",
  varname_date = "TIMESTAMP",
  thresh_drop = 0.0
  )

Retain only events of a minimum length of 20 days.

out_cwd$inst <- out_cwd$inst |>
  filter(len >= 20)

Plot CWD time series.

ggplot() +
  geom_rect(
    data = out_cwd$inst,
    aes(xmin = date_start, xmax = date_end, ymin = -99, ymax = 99999),
    fill = rgb(0,0,0,0.3),
    color = NA) +
  geom_line(data  =  out_cwd$df, aes(TIMESTAMP, prec), size  =  0.3, color = "royalblue") +
  geom_line(data  =  out_cwd$df, aes(TIMESTAMP, deficit), color = "tomato") +
  coord_cartesian(ylim = c(0, 170)) +
  theme_classic() +
  labs(x = "Date", y = "Cumulative water deficit (mm)")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Event analysis

Plot any variable contained in the original data frame as a “superimposed epoch analysis”, that is with data from multiple events (indicated by index iinst) versus the day into event (dday).

out_cwd$df |> 
  ggplot(aes(dday, et, group = iinst)) +
  geom_line(alpha = 0.5) +
  theme_classic()

## Warning: Removed 702 rows containing missing values or values outside the scale range
## (`geom_line()`).

The instances (events) are available in summarised form:

out_cwd$inst

##    idx_start len iinst date_start   date_end max_deficit idx_max_deficit
## 1         22  22     4 2005-01-22 2005-02-12    7.764306              39
## 2         44  33     5 2005-02-13 2005-03-17   22.014083              72
## 3        144  91    13 2005-05-24 2005-08-22   76.584127             198
## 4        236  40    14 2005-08-24 2005-10-02   39.072879             251
## 5        278 137    15 2005-10-05 2006-02-18   37.219533             382
## 6        518  88    24 2006-06-02 2006-08-28  106.801695             573
## 7        646  43    29 2006-10-08 2006-11-19   20.232584             666
## 8        825  43    49 2007-04-05 2007-05-17   46.022031             851
## 9        978  23    65 2007-09-05 2007-09-27   19.154034             990
## 10      1001  45    66 2007-09-28 2007-11-11   31.778082            1019
## 11      1075  27    69 2007-12-11 2008-01-06   10.759122            1098
## 12      1134  23    72 2008-02-08 2008-03-01    5.028310            1151
## 13      1217 205    80 2008-05-01 2008-11-21  123.280649            1318
## 14      1436  63    82 2008-12-06 2009-02-06   15.863117            1449
## 15      1503  20    84 2009-02-11 2009-03-02   11.271995            1517
## 16      1572  23    90 2009-04-21 2009-05-13   16.475135            1577
## 17      1597 286    91 2009-05-16 2010-02-25  153.669072            1741
## 18      1922  30    94 2010-04-06 2010-05-05   39.106875            1945
## 19      1998 186    98 2010-06-21 2010-12-23  117.693244            2028
## 20      2286  65   105 2011-04-05 2011-06-08   37.333746            2322
## 21      2361  25   107 2011-06-19 2011-07-13   39.448022            2378
## 22      2396  43   110 2011-07-24 2011-09-04   34.688899            2436
## 23      2454  20   112 2011-09-20 2011-10-09   30.211974            2469
## 24      2484  48   114 2011-10-20 2011-12-06   22.007549            2526
## 25      2584  31   123 2012-01-28 2012-02-27   11.610848            2609
## 26      2747  92   138 2012-07-09 2012-10-08   71.835714            2790
## 27      2929  24   148 2013-01-07 2013-01-30   11.648926            2950
## 28      2959  31   150 2013-02-06 2013-03-08   18.555859            2985
## 29      3085 148   163 2013-06-12 2013-11-06  127.542176            3171
## 30      3251  25   168 2013-11-25 2013-12-19   11.041014            3271
## 31      3350  25   180 2014-03-04 2014-03-28   12.542856            3367
## 32      3435  64   186 2014-05-28 2014-07-30   86.621857            3474
## 33      3515  37   191 2014-08-16 2014-09-21   25.587160            3524
## 34      3552  26   192 2014-09-22 2014-10-17   15.612616            3569
## 35      3610  26   198 2014-11-19 2014-12-14   12.422168            3631

You can for example plot the distribution of event (instance) length (here in days):

out_cwd$inst |> 
  ggplot(aes(len, after_stat(count))) +
  geom_histogram(color = "black", fill = "grey60") +
  labs(x = "Event length (d)")

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Extreme value statistics

Get annual maxima and fit a general extreme value distribution using the {extRemes} package.

vals <- out_cwd$inst %>%
  group_by(year(date_start)) %>%
  summarise(deficit = max(max_deficit, na.rm = TRUE)) %>%
  pull(deficit)

evd_gev <- extRemes::fevd(x = vals, type = "GEV", method = "MLE", units = "years")
summary(evd_gev)

## 
## extRemes::fevd(x = vals, type = "GEV", method = "MLE", units = "years")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  49.56828 
## 
## 
##  Estimated parameters:
##  location     scale     shape 
## 84.723272 36.697008 -0.412157 
## 
##  Standard Error Estimates:
##   location      scale      shape 
## 13.2609250 10.3279923  0.2916123 
## 
##  Estimated parameter covariance matrix.
##            location      scale       shape
## location 175.852132  -1.605574 -1.80779489
## scale     -1.605574 106.667425 -2.10522917
## shape     -1.807795  -2.105229  0.08503775
## 
##  AIC = 105.1366 
## 
##  BIC = 106.0443

Get CWD magnitudes for given 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(
  evd_gev,
  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(return_period, return_level)) +
  geom_point() +
  labs(x = "Return period (yr)", 
       y = "Magnitude of annual CWD maximum (mm)",
       title = "GEV")

With a Gumbel extreme value distribution, the return period as a function of the CWD extreme magnitude is calculated as follows:

# Fit Gumbel distribution
evd_gumbi <- extRemes::fevd(x = vals, type = "Gumbel", method = "MLE", units = "years")
summary(evd_gumbi)

## 
## extRemes::fevd(x = vals, type = "Gumbel", method = "MLE", units = "years")
## 
## [1] "Estimation Method used: MLE"
## 
## 
##  Negative Log-Likelihood Value:  50.26535 
## 
## 
##  Estimated parameters:
## location    scale 
## 77.33370 32.62668 
## 
##  Standard Error Estimates:
##  location     scale 
## 10.936112  7.919814 
## 
##  Estimated parameter covariance matrix.
##           location    scale
## location 119.59854 28.88856
## scale     28.88856 62.72345
## 
##  AIC = 104.5307 
## 
##  BIC = 105.1359

# calculate return period as a function of the CWD extreme. Using the two 
# coefficients of the fitted distribution as arguments
calc_return_period <- function(x, loc, scale){
  1 / (1 - exp(-exp(-(x-loc)/scale)))
}

extract_loc <- function(mod){
  loc <- mod$results$par[ "location" ]
  if (!is.null(loc)){
    return(loc)
  } else {
    return(NA)
  }
}

extract_scale <- function(mod){
  scale <- mod$results$par[ "scale" ]
  if (!is.null(scale)){
    return(scale)
  } else {
    return(NA)
  }
}

# demo return periods
return_period <- c(2, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 200, 250, 300, 500, 800)

# use built-in function to get expected CWD extreme for given return periods 
# (inverse of probability)
return_level <- extRemes::return.level(
  evd_gumbi,
  return.period = return_period
)

# create data frame for visualisation
df_return <- tibble(
  return_period = return_period,
  return_level = unname(c(return_level)),
  trans_level = -log( -log(1 - 1/return_period))) |> 
  mutate(myreturn_period = calc_return_period(
    return_level,
    extract_loc(evd_gumbi),
    extract_scale(evd_gumbi)
  ))

# CWD extreme for a given return period
df_return |> 
  ggplot(aes(return_period, return_level)) +
  geom_point() +
  labs(x = "Return period (yr)", 
       y = "Magnitude of annual CWD maximum (mm)",
       title = "Gumbel")

# Return period for a given CWD extreme (calculated based on function above)
df_return |> 
  ggplot(aes(return_level, myreturn_period)) + 
  geom_point() +
  labs(y = "Return period (yr)", 
       x = "Magnitude of annual CWD maximum (mm)",
       title = "Gumbel")

Visualise the estimated event size with a return period of $T = 80$ y as the red line on top of the distribution of cumulative water deficit events.

ggplot() +
  geom_histogram(
    data = out_cwd$inst,
    aes(x = max_deficit, y = after_stat(density)),
    color = "black",
    position="identity",
    bins = 6
    ) +
  labs(x = "Cumulative water deficit (mm)") +
  geom_vline(xintercept = df_return %>%
               dplyr::filter(return_period == 80) %>%
               pull(return_level),
             col = "tomato")

Time stepping

The data frame used above contains time series with daily resolution. The CWD algorithm can also be applied to data provided at other time steps. It primarily acts on the rows in the data frame.

wdf <- df |> 
  mutate(year = lubridate::year(TIMESTAMP),
         week = lubridate::week(TIMESTAMP)) |> 
  group_by(year, week) |> 
  summarise(wbal = sum(wbal, na.rm = FALSE)) |> 
  
  # create a date object again, considering the first day of the week
  mutate(date = lubridate::ymd(paste0(year, "-01-01")) + lubridate::weeks(week-1))

## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by year and week.
## ℹ Output is grouped by year.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(year, week))` for per-operation grouping
##   (`?dplyr::dplyr_by`) instead.

out_cwd_weekly <- cwd(
  wdf,
  varname_wbal = "wbal",
  varname_date = "date",
  thresh_drop = 0.0
  )

Plot weekly CWD time series in green (Daily CWD time series are plotted by the red line.)

ggplot() +
  geom_rect(
    data = out_cwd_weekly$inst,
    aes(xmin = date_start, xmax = date_end, ymin = -99, ymax = 99999),
    fill = rgb(0,0,0,0.3),
    color = NA) +
  geom_line(data  =  out_cwd$df, aes(TIMESTAMP, deficit), color = "tomato", alpha = 0.6) +
  geom_line(data  =  out_cwd_weekly$df, aes(date, deficit), color = "springgreen4") +
  coord_cartesian(ylim = c(0, 170)) +
  theme_classic() +
  labs(x = "Date", y = "Cumulative water deficit (mm)")

Beni Stocker

2026-07-29