This markdown script reproduces the supplementary figures associated with the paper “Atmospheric River Sequences as Indicators of Hydrologic Hazard in Historical Reanalysis and GFDL SPEAR Future Climate Projections” (https://doi.org/10.22541/essoar.167590838.86645781/v1).

Setup

Load packages & functions

## set random seed to produce consistent results
set.seed(2023)

source('_data/setup.R')
source('_scripts/create_df_functions.R')

## additional packages not in setup.R
library(patchwork)
library(pammtools)
library(lhs)

## progress bar
progress <- FALSE

## geospatial helpers
colors.huc <- colors.huc[c(6:10,1:5)]

## download custom boxplot function
library(assertthat)
library(devtools)
source_url('https://github.com/BoulderCodeHub/CRSSIO/blob/master/R/stat_boxplot_custom.R?raw=TRUE')

Load MERRA-2 catalogs

# ## MERRA 3-hour dataset
# load('_scripts/_checkpoints/df_3hr_1209.Rdata')

## MERRA 24-hour dataset
load('_scripts/_checkpoints/df_24hr_0209.Rdata')

## MERRA metadata 
ts_merra <- seq(ymd_hms('1980-1-1 00:00:00'), ymd_hms('2021-12-31 21:00:00'), by = '3 hours')
dates_merra <- unique(as.Date(ts_merra))

Load GFDL catalogs

## SPEAR historic
load('_scripts/_checkpoints/df_hist_1212.Rdata')

## SPEAR SSP 2-4.5
load('_scripts/_checkpoints/df_ssp245_1214.Rdata')

## SPEAR SSP 5-8.5
load('_scripts/_checkpoints/df_ssp585_1214.Rdata')

## SPEAR metadata
load('_data/GFDL/gfdl_metadata.Rdata')  #ts_hist, ts_future
ts_decades <- data.frame(ts = ts_future) %>% 
  filter(year(ts) > 2020 & year(ts) <= 2090) %>% 
  pull(ts)
decades <- 
  data.frame(yr = 2021:2090) %>% 
  mutate(decade = ceiling(yr/10)-202)

FIG S1: Validation of bias correction method

compare_biascorr <-
  map_dfr(.x = 1801:1810, .f = function(huc) {
    i <- which(grid_hucrep[] == huc)
    df_hist[[i]] %>%
      filter(wateryear(date) %in% 1981:2010) %>%
      select(date,ens,ivt) %>%
      pivot_wider(names_from = ens, names_prefix = 'ens0', values_from = ivt) %>%
      select(-date) %>%
      mutate(
        merra = df_24hr[[i]] %>%
          filter(wateryear(date) %in% 1981:2010) %>%
          pull(ivt.snap) %>% sort) %>%
      apply(2, sort) %>% as.data.frame %>%
      pivot_longer(starts_with('ens'), names_to = 'ens', values_to = 'gfdl') %>%
      mutate(
        huc4 = paste0('HUC', huc), 
        tag = c(paste0('(',letters[huc-1800],')'), rep(NA, nrow(.)-1)))
  }) %>% filter(merra < 2000 & gfdl < 2000)

ggplot(compare_biascorr) + 
  geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) + 
  geom_step(aes(x = merra, y = gfdl, group = ens, color = factor(ens)), direction = 'mid') + 
  geom_text(
    aes(x = 250, y = 1200, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') + 
  facet_wrap(~huc4, nrow = 2) + 
  scale_color_manual('Ensemble Member', labels = 1:5, values = paste0('grey', seq(20,80,15))) + 
  guides(color = guide_legend(override.aes = list(size = 1))) + 
  scale_x_origin('MERRA-2 IVT (kg/m/s)', breaks = 250*(0:10)) + 
  scale_y_origin('GFDL SPEAR Bias-Corrected IVT (kg/m/s)', breaks = 250*(0:10)) + 
  geom_parity() + coord_fixed() +
  theme(
    text = element_text(size = 10),
    panel.grid.major = element_line(size = 0.25),
    axis.line = element_blank(),
    legend.position = 'bottom',
    strip.background = element_rect(fill = 'grey95', size = NA),
    strip.text = element_text(margin = margin(2,2,2,2)),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s1_biascorr.png', width = 6, height = 3.75, dpi = 600)

FIG S2: Comparison of MERRA-2 versus bias-corrected GFDL SPEAR annual average sequence day frequency

compare_freq <-
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i) 
    if (i %in% index_ca) {
      temp <- df_hist[[i]] %>% 
        filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)) %>% 
        count(ens,seq) %>% filter(seq)
      merra <- df_24hr[[i]] %>% 
        filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)) %>% 
        count(seq) %>% filter(seq) %>% pull(n)
      temp %>% mutate(merra = if (length(merra)==0) NA else merra) %>% mutate(id = i)
    }
  }

r2 <- map(
    .x = 1:5,
    .f = ~summary(lm(merra ~ n, data = compare_freq %>% filter(ens==.x)))$r.squared) %>% 
  reduce(c)

compare_freq %>% 
  group_by(ens) %>% 
  mutate(
    id = 1:length(id),
    tag = case_when(id==1 ~ paste0('(', letters[ens], ')')),
    rsq = case_when(id==1 ~ paste0('R^2','==', round(r2,3)[ens]))) %>% 
  ggplot() + 
  geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) + 
  geom_point(aes(x = merra/30, y = n/30), size = 1) + 
  geom_text(
    aes(x = 10, y = 81, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') + 
  geom_text(
    aes(x = 80, y = 20, label = rsq), hjust = 1, parse = TRUE,
    family = 'Segoe UI', size = 9/.pt) + 
  facet_wrap(~paste('Ensemble Member',ens)) + 
  scale_x_origin(
    'MERRA-2 Annual Average Sequence Days', breaks = 25*(0:5)) + 
  scale_y_origin(
    'GFDL SPEAR Annual Average Sequence Days', breaks = 25*(0:5)) + 
  geom_abline(linetype = 'dashed', color = 'gray30') + coord_fixed() + 
  theme(
    text = element_text(size = 10),
    panel.grid.major = element_line(size = 0.25),
    legend.position = 'bottom',
    axis.line = element_blank(),
    strip.background = element_rect(fill = 'grey95', size = NA),
    strip.text = element_text(margin = margin(2,2,2,2)),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
 ggsave('_figures/s2_biascorr_freq.png', width = 6, height = 4.75, dpi = 600)

FIG S3: Comparison of MERRA-2 versus bias-corrected GFDL SPEAR sequence maximum IVT

compare <-
  foreach (huc = 1801:1810) %do% {
    i <- which(grid_hucrep[] == huc)
    rbind(
      df_hist[[i]] %>% 
        filter(seq) %>% group_by(seq.count) %>% 
        summarize(date = date[1], maxivt = seq.maxivt[1], duration = seq.duration[1], ens = ens[1]) %>% 
        filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)),
      df_24hr[[i]] %>% 
        filter(seq) %>% group_by(seq.count) %>% 
        summarize(date = date[1], maxivt = seq.maxivt[1], duration = seq.duration[1], ens = NA) %>% 
        filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3))) %>% 
      mutate(huc = huc)
  } %>% 
  reduce(rbind) %>% 
  group_by(huc) %>% 
  mutate(
    id = 1:length(huc),
    tag = case_when(id==1 ~ paste0('(',letters[huc-1800],')'))) %>% 
  ungroup
compare %>% 
  group_by(ens,huc) %>% 
  arrange(maxivt) %>% 
  mutate(p = add_index(length(maxivt))) %>% 
  ungroup %>% 
  ggplot() + 
  geom_vline(xintercept = 250, size = 0.75) + geom_hline(yintercept = 0) + 
  geom_step(
    aes(x = maxivt, y = p, group = ens, color = is.na(ens), size = is.na(ens))) + 
  geom_text(
    aes(x = 380, y = 0.9, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') + 
  facet_wrap(~paste0('HUC',huc), nrow = 2) +
  scale_color_manual(
    'Data Source',
    values = c('grey60', 'black'),
    labels = c('GFDL SPEAR', 'MERRA-2')) + 
  scale_size_manual(values = c(0.25, 0.5)) + 
  guides(
    color = guide_legend(override.aes = list(size = 1)),
    size = guide_none()) + 
  scale_x_continuous(
    'Sequence Max IVT (kg/m/s)', breaks = 250*(0:10),
    expand = expansion(mult = c(0,0.05))) + 
  scale_y_origin('Cumulative Probability') +
  theme(
    text = element_text(size = 10),
    panel.grid.major = element_line(size = 0.25),
    legend.position = 'bottom',
    strip.background = element_rect(fill = 'grey95', size = NA),
    strip.text = element_text(margin = margin(2,2,2,2)),
    axis.line = element_blank(),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s3_biascorr_maxivt.png', width = 6, height = 3.75, dpi = 600)

FIG S4: Comparison of MERRA-2 versus bias-corrected GFDL SPEAR sequence duration

compare %>% 
  group_by(ens,huc) %>% 
  arrange(duration) %>% 
  mutate(p = add_index(length(duration))) %>% 
  ungroup %>% 
  ggplot() + 
  geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) + 
  geom_step(
    aes(x = duration, y = p, group = ens, color = is.na(ens), size = is.na(ens))) + 
  geom_text(
    aes(x = 13, y = 0.9, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') + 
  facet_wrap(~paste0('HUC',huc), nrow = 2) +
  scale_color_manual(
    'Data Source',
    values = c('grey60', 'black'),
    labels = c('GFDL SPEAR', 'MERRA-2')) + 
  scale_size_manual(values = c(0.25, 0.5)) + 
  guides(
    color = guide_legend(override.aes = list(size = 1)),
    size = guide_none()) + 
  scale_x_origin('Sequence Duration (days)', breaks = 30*(0:10)) + 
  scale_y_origin('Cumulative Probability') +
  theme(
    text = element_text(size = 10),
    panel.grid.major = element_line(size = 0.25),
    legend.position = 'bottom',
    strip.background = element_rect(fill = 'grey95', size = NA),
    strip.text = element_text(margin = margin(2,2,2,2)),
    axis.line = element_blank(),
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s4_biascorr_duration.png', width = 6, height = 3.75, dpi = 600)

FIG S5: Sensitivity analysis results

Create clustered parameters & outcomes

Generate LHS samples

## define big table of samples
samples.norm <- improvedLHS(n = 100, k = 3)
samples.dist <- 
  data.frame(
    window = samples.norm[,1] %>% qunif(min = 3, max = 10) %>% round,
    startend = samples.norm[,2] %>% qunif(min = 0.25, max = 0.75),
    magcutoff = samples.norm[,3] %>% qunif(min = 200, max = 300)) %>% 
  mutate(magcutoff = cbind(magcutoff, startend) %>% apply(1, max))
samples.dist <- samples.dist %>% rbind(c(window = 5, startend = 0.5, magcutoff = 250))

Generate output statistics for every LHS parameter combination

(this takes a while \(\to\) load from file)

# ## define progress bar
# if (progress) pb <- txtProgressBar(min = 0, max = 10*nrow(samples.dist), style = 3)

# ## historical
# timer <- Sys.time()
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_hist <- 
#   foreach(
#     huc = 1801:1810,
#     .options.snow = if (progress) opts,
#     .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
#     .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
#   foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% { 
#     ## create sequence catalog
#     data <- df_hist[[which(grid_hucrep[]==huc)]] %>% 
#       filter(wateryear(date) %in% 1981:2010) %>% 
#       transmute(date, ivt, ens)
#     sequences <-
#       map(.x = 1:5, .f = function(x) {
#         if (nrow(data %>% filter(ens==x)) > 0) {
#           calculate_sequences(
#             data %>% filter(ens==x), data %>% filter(ens==x), 
#             window = samples.dist$window[samp], 
#             startend = samples.dist$startend[samp], 
#             magcutoff = samples.dist$magcutoff[samp]) %>% 
#             mutate(ens=x)
#         }
#       }) %>% reduce(rbind)
#     ## calculate sequence statistics
#     sequences %>%
#       group_by(ens) %>% 
#         mutate(
#           nseq = length(ens),
#           freq = nseq/30) %>% 
#         ungroup %>% 
#         summarize(
#           maxivt.mean = mean(maxivt),
#           maxivt.q95 = quantile(maxivt,0.95),
#           duration.mean = mean(duration),
#           duration.q95 = quantile(duration,0.95),
#           nseq = nseq[1], freq = freq[1])
#   }
# stopCluster(clust)
# if (progress) cat('\n')

# ## SSP 2-4.5
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_ssp245 <- 
#   foreach(
#     huc = 1801:1810,
#     .options.snow = if (progress) opts,
#     .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
#     .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
#   foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% { 
#     ## create sequence catalog
#     data <- df_ssp245[[which(grid_hucrep[]==huc)]] %>% 
#       filter(wateryear(date) %in% 2061:2090) %>% 
#       transmute(date, ivt, ens)
#     data_hist <- df_hist[[which(grid_hucrep[]==huc)]] %>% 
#       filter(wateryear(date) %in% 1981:2010) %>% 
#       transmute(date, ivt, ens)
#     sequences <-
#       map(.x = 1:5, .f = function(x) {
#         if (nrow(data %>% filter(ens==x)) > 0) {
#           calculate_sequences(
#             data %>% filter(ens==x), data_hist %>% filter(ens==x), 
#             window = samples.dist$window[samp], 
#             startend = samples.dist$startend[samp], 
#             magcutoff = samples.dist$magcutoff[samp]) %>% 
#             mutate(ens=x)
#         }
#       }) %>% reduce(rbind)
#     ## calculate sequence statistics
#     sequences %>%
#       group_by(ens) %>% 
#         mutate(
#           nseq = length(ens),
#           freq = nseq/30) %>% 
#         ungroup %>% 
#         summarize(
#           maxivt.mean = mean(maxivt),
#           maxivt.q95 = quantile(maxivt,0.95),
#           duration.mean = mean(duration),
#           duration.q95 = quantile(duration,0.95),
#           nseq = nseq[1], freq = freq[1])
#   }
# stopCluster(clust)
# if (progress) cat('\n')

# ## SSP 5-8.5
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_ssp585 <- 
#   foreach(
#     huc = 1801:1810,
#     .options.snow = if (progress) opts,
#     .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
#     .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
#   foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% { 
#     ## create sequence catalog
#     data <- df_ssp585[[which(grid_hucrep[]==huc)]] %>% 
#       filter(wateryear(date) %in% 2061:2090) %>% 
#       transmute(date, ivt, ens)
#     data_hist <- df_hist[[which(grid_hucrep[]==huc)]] %>% 
#       filter(wateryear(date) %in% 1981:2010) %>% 
#       transmute(date, ivt, ens)
#     sequences <-
#       map(.x = 1:5, .f = function(x) {
#         if (nrow(data %>% filter(ens==x)) > 0) {
#           calculate_sequences(
#             data %>% filter(ens==x), data_hist %>% filter(ens==x), 
#             window = samples.dist$window[samp], 
#             startend = samples.dist$startend[samp], 
#             magcutoff = samples.dist$magcutoff[samp]) %>% 
#             mutate(ens=x)
#         }
#       }) %>% reduce(rbind)
#     ## calculate sequence statistics
#     sequences %>%
#       group_by(ens) %>% 
#         mutate(
#           nseq = length(ens),
#           freq = nseq/30) %>% 
#         ungroup %>% 
#         summarize(
#           maxivt.mean = mean(maxivt),
#           maxivt.q95 = quantile(maxivt,0.95),
#           duration.mean = mean(duration),
#           duration.q95 = quantile(duration,0.95),
#           nseq = nseq[1], freq = freq[1])
#   }
# stopCluster(clust)
# if (progress) cat('\n')
# Sys.time() - timer

# ## save checkpoint
# save(samples.dist, sens_hist, sens_ssp245, sens_ssp585, file = '_results/sensitivity_checkpoint.Rdata')
## load checkpoint from file
load('_results/_checkpoints/sensitivity_checkpoint.Rdata')

Combine outcomes into dataframes

sens_hist <- 
  lapply(1:10, function(i) sens_hist[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>% 
  reduce(rbind) %>% 
  mutate(freq = setNA(freq,0))
sens_ssp245 <- 
  lapply(1:10, function(i) sens_ssp245[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>% 
  reduce(rbind) %>% 
  mutate(freq = setNA(freq,0))
sens_ssp585 <- 
  lapply(1:10, function(i) sens_ssp585[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>% 
  reduce(rbind) %>% 
  mutate(freq = setNA(freq,0))

## outcome: absolute value
sens <- 
  rbind(
    sens_hist %>% mutate(source = 'hist'),
    sens_ssp245 %>% mutate(source = 'ssp245'),
    sens_ssp585 %>% mutate(source = 'ssp585')) %>% 
  left_join(samples.dist %>% mutate(index = 1:nrow(.)), by = 'index')

## outcome: relative change
sens_change_relative <- 
  rbind(
    (sens_ssp245 / sens_hist) %>% select(-huc, -index) %>% 
      cbind(sens_hist %>% select(huc, index)) %>% 
      mutate(source = 'ssp245'),
    (sens_ssp585 / sens_hist) %>% select(-huc, -index) %>% 
      cbind(sens_hist %>% select(huc, index)) %>% 
      mutate(source = 'ssp585')) %>% 
  left_join(samples.dist %>% mutate(index = 1:nrow(.)), by = 'index')

Create clusters

## decide on an appropriate number of clusters
k <- 3

## define custom kmeans function
Kmeans <- function(x, ...) kmeans(x, centers = k, nstart = 50, ...)

## outcome: absolute value
sens <- 
  foreach (h = 1801:1810, .combine = 'rbind') %do% { 
    temp <- sens %>% filter(huc == h)
    good.freq <- 
      !is.na(temp$freq) & !is.nan(temp$freq) & !is.infinite(temp$freq)
    good.maxivt <- 
      !is.na(temp$maxivt.mean) & !is.nan(temp$maxivt.mean) & !is.infinite(temp$maxivt.mean)
    good.duration <- 
      !is.na(temp$duration.mean) & !is.nan(temp$duration.mean) & !is.infinite(temp$duration.mean)
    temp %>% 
      mutate(
        cl.freq = ifelse(
          good.freq, Kmeans(freq[good.freq])$cluster, NA),
        cl.maxivt = ifelse(
          good.maxivt, Kmeans(maxivt.mean[good.maxivt])$cluster, NA),
        cl.duration = ifelse(
          good.duration, Kmeans(duration.mean[good.duration])$cluster, NA))
  }

## outcome: relative change
sens_change_relative <- 
  foreach (h = 1801:1810, .combine = 'rbind') %do% { 
    temp <- sens_change_relative %>% filter(huc == h)
    good.freq <- 
      !is.na(temp$freq) & !is.nan(temp$freq) & !is.infinite(temp$freq)
    good.maxivt <- 
      !is.na(temp$maxivt.mean) & !is.nan(temp$maxivt.mean) & !is.infinite(temp$maxivt.mean)
    good.duration <- 
      !is.na(temp$duration.mean) & !is.nan(temp$duration.mean) & !is.infinite(temp$duration.mean)
    temp %>% 
      mutate(
        cl.freq = ifelse(
          good.freq, Kmeans(freq[good.freq])$cluster, NA),
        cl.maxivt = ifelse(
          good.maxivt, Kmeans(maxivt.mean[good.maxivt])$cluster, NA),
        cl.duration = ifelse(
          good.duration, Kmeans(duration.mean[good.duration])$cluster, NA))
  }

Perform response sensitivity analysis (RSA)

Define L1 bootstrap function

## bootstrap L1-norm distribution
L1_boot <- function(data, cl, var, boot = 1e3) {
  ## find existing parameters within clusters
  x <- sort(data[data$cl == cl,] %>% pull(get(var)))
  x.all <- sort(data %>% pull(get(var)))
  p <- (1:length(x))/length(x)
  p.all <- (1:length(x.all))/length(x.all)
  x.interp <- interp1(x = p.all, y = x.all, xi = p)

  if (is.na(boot)) {
    return(data.frame(x, x.interp) %>% as.matrix %>% apply(1, diff) %>% abs %>% sum)
  } else {
    L1 <- c()
    for (b in 1:boot) {
      L1[b] <- 
        data.frame(
          x.boot = sort(sample(x.all, size = length(x), replace = TRUE)), 
          x.interp) %>% 
        as.matrix %>% apply(1, diff) %>% abs %>% sum
    }
    return(L1)
  }
}

Calculate bootstrapped likelihood of L1 distance


Frequency:

if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
clust <- parallel::makeCluster(cores)
registerDoSNOW(clust)
rsa_freq <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens %>%
        filter(huc==h & !is.na(cl.freq)) %>% rename(cl = cl.freq)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>%
        cbind(
          L1_boot(temp, cl, var) %>%
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
cat('\n')

rsa_freq_rel <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens_change_relative %>%
        filter(huc==h & !is.na(cl.freq)) %>% rename(cl = cl.freq)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>%
        cbind(
          L1_boot(temp, cl, var) %>%
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)


Intensity:

if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
clust <- parallel::makeCluster(cores)
registerDoSNOW(clust)
rsa_maxivt <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens %>% 
        filter(huc==h & !is.na(cl.maxivt)) %>% rename(cl = cl.maxivt)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>% 
        cbind(
          L1_boot(temp, cl, var) %>% 
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
cat('\n')

rsa_maxivt_rel <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens_change_relative %>% 
        filter(huc==h & !is.na(cl.maxivt)) %>% rename(cl = cl.maxivt)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>% 
        cbind(
          L1_boot(temp, cl, var) %>% 
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)


Duration:

if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
clust <- parallel::makeCluster(cores)
registerDoSNOW(clust)
rsa_duration <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens %>% 
        filter(huc==h & !is.na(cl.duration)) %>% rename(cl = cl.duration)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>% 
        cbind(
          L1_boot(temp, cl, var) %>% 
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
cat('\n')

rsa_duration_rel <- 
  foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:% 
  foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
  foreach (
    var = c('startend', 'window', 'magcutoff'), 
    .options.snow = if (progress) opts,
    .export = 'L1_boot',
    .packages = c('tidyverse', 'pracma'),
    .combine = 'rbind', .inorder = FALSE) %dopar% {
      print(c(h, cl, var))
      ## generate normalized L1 distances for every parameter & cluster
      temp <- sens_change_relative %>% 
        filter(huc==h & !is.na(cl.duration)) %>% rename(cl = cl.duration)
      data.frame(
        cl, var, huc = h,
        d = L1_boot(temp, cl, var, NA)) %>% 
        cbind(
          L1_boot(temp, cl, var) %>% 
            quantile(c(0.95, 0.99, 0.995)) %>%
            t %>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
    } %>% 
  select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)

Create sensitivity analysis plot

Create plot panels by row


Frequency:

g1.freq <- ggplot(sens) + 
  geom_smooth(
    aes(x = jitter(window), y = freq, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.x = element_blank())
g4.freq <- ggplot(sens) + 
  geom_smooth(
    aes(x = startend, y = freq, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())
g7.freq <- ggplot(sens) + 
  geom_smooth(
    aes(x = magcutoff, y = freq, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())

lims <- 
  rbind(
    layer_scales(g1.freq)$y$get_limits(),
    layer_scales(g4.freq)$y$get_limits(),
    layer_scales(g7.freq)$y$get_limits())
lims <- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
g1.freq <- g1.freq + coord_cartesian(ylim = lims)
g4.freq <- g4.freq + coord_cartesian(ylim = lims)
g7.freq <- g7.freq + coord_cartesian(ylim = lims)

g1.rsa <- rsa_freq %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  
g3.freq <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = jitter(window), y = freq-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  scale_y_continuous(labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.x = element_blank())
g6.freq <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = startend, y = freq-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())
g9.freq <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = magcutoff, y = freq-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())

lims <- 
  rbind(
    layer_scales(g3.freq)$y$get_limits(),
    layer_scales(g6.freq)$y$get_limits(),
    layer_scales(g9.freq)$y$get_limits())
lims <- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
g3.freq <- g3.freq + coord_cartesian(ylim = lims)
g6.freq <- g6.freq + coord_cartesian(ylim = lims)
g9.freq <- g9.freq + coord_cartesian(ylim = lims)

g2.rsa <- rsa_freq_rel %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  


Intensity:

g1.maxivt <- ggplot(sens) + 
  geom_smooth(
    aes(x = jitter(window), y = maxivt.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.x = element_blank())
g4.maxivt <- ggplot(sens) + 
  geom_smooth(
    aes(x = startend, y = maxivt.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())
g7.maxivt <- ggplot(sens) + 
  geom_smooth(
    aes(x = magcutoff, y = maxivt.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())

lims <- 
  rbind(
    layer_scales(g1.maxivt)$y$get_limits(),
    layer_scales(g4.maxivt)$y$get_limits(),
    layer_scales(g7.maxivt)$y$get_limits())
lims <- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
g1.maxivt <- g1.maxivt + coord_cartesian(ylim = lims)
g4.maxivt <- g4.maxivt + coord_cartesian(ylim = lims)
g7.maxivt <- g7.maxivt + coord_cartesian(ylim = lims)

g3.rsa <- rsa_maxivt %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  
g3.maxivt <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = jitter(window), y = maxivt.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.x = element_blank())
g6.maxivt <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = startend, y = maxivt.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())
g9.maxivt <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = magcutoff, y = maxivt.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())

g4.rsa <- rsa_maxivt_rel %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity', breaks = 0:3) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  


Duration:

g1.duration <- ggplot(sens) + 
  geom_smooth(
    aes(x = jitter(window), y = duration.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.x = element_blank())
g4.duration <- ggplot(sens) + 
  geom_smooth(
    aes(x = startend, y = duration.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())
g7.duration <- ggplot(sens) + 
  geom_smooth(
    aes(x = magcutoff, y = duration.mean, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text = element_blank())

lims <- 
  rbind(
    layer_scales(g1.duration)$y$get_limits(),
    layer_scales(g4.duration)$y$get_limits(),
    layer_scales(g7.duration)$y$get_limits())
lims <- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
g1.duration <- g1.duration + coord_cartesian(ylim = lims)
g4.duration <- g4.duration + coord_cartesian(ylim = lims)
g7.duration <- g7.duration + coord_cartesian(ylim = lims)

g5.rsa <- rsa_duration %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  
g3.duration <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = jitter(window), y = duration.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(text = element_text(family = 'Segoe UI', size = 8), axis.title = element_blank())
g6.duration <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = startend, y = duration.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.y = element_blank())
g9.duration <- ggplot(sens_change_relative) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.25) + 
  geom_smooth(
    aes(x = magcutoff, y = duration.mean-1, group = huc, color = factor(huc)), 
    size = 0.75, se = FALSE) + 
  scale_color_manual('', values = colors.huc, guide = guide_none()) + 
  scale_y_continuous(labels = percent, limits = c(0,1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    axis.title = element_blank(), axis.text.y = element_blank())

g6.rsa <- rsa_duration_rel %>% 
  mutate(
    var_name = case_when(
      var == 'window' ~ 'Moving\nAverage',
      var == 'startend' ~ 'Start-End\nThreshold',
      var == 'magcutoff' ~ 'Mag.\nFilter') %>% 
      factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>% 
  group_by(var_name, huc) %>% 
  summarize(d.norm = mean(d/d95), .groups = 'drop') %>% 
  ggplot() + 
  geom_col(
    aes(
      y = d.norm, x = var_name, 
      color = d.norm>1, alpha = d.norm>1, 
      group = factor(huc), fill = factor(huc)),
    position = 'dodge', size = 0.25) + 
  geom_hline(yintercept = 1, linetype = 'dashed') + 
  scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) + 
  scale_fill_manual('', values = colors.huc, guide = guide_none()) + 
  scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) + 
  scale_y_origin('Sensitivity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title.x = element_blank(), legend.position = 'bottom')  

Combine plot panels into a single figure

Note: Additional labels were added to this figure in Inkscape.

ggplot(sens) + 
  geom_histogram(aes(x = huc, fill = factor(huc))) + 
  scale_fill_manual(
    'Hydrologic Region', 
    values = colors.huc, guide = guide_legend(direction = 'horizontal', nrow = 1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8), 
    legend.position = 'bottom',
    legend.text = element_text(margin = margin(-2,-2,-2,-2)),
    legend.margin = margin(-2,-2,-2,-2), 
    legend.box.margin = margin(-2,-2,-2,-2))
ggsave('_figures/s5_legend.png', width = 5.5, height = 1, dpi = 600)

g1.freq + g4.freq + g7.freq + g1.rsa + theme(axis.text.x = element_blank()) + 
  g3.freq + g6.freq + g9.freq + g2.rsa + theme(axis.text.x = element_blank()) + 
  g1.maxivt + g4.maxivt + g7.maxivt + g3.rsa + theme(axis.text.x = element_blank()) + 
  g3.maxivt + g6.maxivt + g9.maxivt + g4.rsa + theme(axis.text.x = element_blank()) + 
  g1.duration + g4.duration + g7.duration + g5.rsa + theme(axis.text.x = element_blank()) + 
  g3.duration + g6.duration + g9.duration + g6.rsa + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  plot_layout(design = 'abcd\nefgh\nijkl\nmnop\nqrst\nuvwx') + 
  plot_annotation(tag_levels = 'a', tag_prefix = '(', tag_suffix = ')') &
  theme(plot.tag = element_text(
    family = 'Segoe UI', size = 10, face = 'bold', margin = margin(-2,0,-2,-10)))
ggsave('_figures/s5.png', width = 5.5, height = 7.75, dpi = 600)

FIG S6: Additive decadal change in number of sequence days

ndays_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (i %in% index_ca) {
      df_hist[[i]] %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        filter(month(date) %in% c(10:12,1:3)) %>% 
        filter(seq) %>% 
        count(ens) %>% 
        summarize(nn = sum(n)/length(unique(ens))/30) %>% 
        pull(nn)
    } else NA
  }
cat('\n')

ndays_ssp245 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp245[[i]])) rep(NA,7) else {
      df_ssp245[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(month(date) %in% c(10:12,1:3)) %>% 
        filter(seq) %>% 
        count(ens, decade) %>% 
        group_by(decade) %>% 
        summarize(nn = sum(n)/length(unique(ens))/10) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        mutate(nn = setNA(nn,0)) %>% 
        pull(nn)
    }
  }
cat('\n')

ndays_ssp585 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp585[[i]])) rep(NA,7) else {
      df_ssp585[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(month(date) %in% c(10:12,1:3)) %>% 
        filter(seq) %>% 
        count(ens, decade) %>% 
        group_by(decade) %>% 
        summarize(nn = sum(n)/length(unique(ens))/10) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        mutate(nn = setNA(nn,0)) %>% 
        pull(nn)
    }
  }
cat('\n')

ndays <- 
  left_join(
    ndays_ssp245 %>% as.data.frame %>%
      cbind(hist = ndays_ref) %>% 
      cbind(huc = grid_hucrep[]) %>% 
      filter(huc > 0) %>% 
      pivot_longer(cols = starts_with('V'), names_to = 'decade') %>% 
      rename(ssp245 = value),
    ndays_ssp585 %>% as.data.frame %>%
      cbind(huc = grid_hucrep[]) %>% 
      filter(huc > 0) %>% 
      pivot_longer(cols = starts_with('V'), names_to = 'decade') %>% 
      rename(ssp585 = value),
    by = c('huc', 'decade')) %>% 
  mutate(
    decade = toNumber(gsub('V', '', decade)),
    decade_name = case_when(
      decade==1 ~ '2020s',
      decade==2 ~ '2030s',
      decade==3 ~ '2040s',
      decade==4 ~ '2050s',
      decade==5 ~ '2060s',
      decade==6 ~ '2070s',
      decade==7 ~ '2080s') %>% factor) %>% 
  pivot_longer(cols = c(ssp245, ssp585), names_to = 'sim') %>% 
  mutate(
    sim_name = case_when(
      sim == 'ssp245' ~ 'SSP 2-4.5 ("Middle of the Road")',
      sim == 'ssp585' ~ 'SSP 5-8.5 ("Upper End")') %>% 
      factor %>% factor(levels = levels(.))) %>%
  mutate(
    id = 1:nrow(.),
    tag = ifelse(id %in% 1:2, paste0('(', letters[id], ')'), NA))
gdays <- ggplot(ndays) + 
  geom_hline(yintercept = 0, color = 'grey50', size = 0.35) +
  geom_line(
    aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)), 
    size = 0.75) + 
  geom_point(
    aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)), 
    size = 1) + 
  geom_text(
    aes(x = 1.5, y = 40, label = tag),
    size = 10/.pt, family = 'Segoe UI', fontface = 'bold') +
  facet_wrap(~sim_name, nrow = 2) + 
  scale_color_manual(values = colors.huc, guide = guide_none()) + 
  scale_x_discrete(expand = expansion(mult = 0.035)) + 
  scale_y_continuous(
    'Change in Number of Sequence Days\nRelative to WY 1981-2010',
    expand = expansion(mult = 0.1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 10),
    panel.grid.major.y = element_line(size = 0.25),
    strip.background = element_rect(color = NA, fill = 'grey95'),
    strip.text = element_text(size = 10, margin = margin(2,0,2,0)),
    legend.position = 'bottom', 
    axis.line.x = element_blank(),
    axis.title.x = element_blank())

gmap <- ggplot(huc4) + 
  geom_sf(aes(fill = factor(huc4)), alpha = 0.5, size = 0.25) + 
  scale_fill_manual('Hydrologic Region', values = colors.huc) + 
  geom_tile(
    data = raster.df(grid_ca) %>% filter(!is.na(value)), 
    aes(x=x, y=y), color = 'grey70', fill = NA) + 
  geom_tile(
    data = raster.df(grid_hucrep) %>% filter(value>0),
    aes(x=x, y=y, fill = factor(value)), color = 'black') + 
  geom_point(
    data = raster.df(grid_hucrep) %>% filter(value>0),
    aes(x=x, y=y), color = 'black', size = 0.5) +
  guides(fill = guide_legend(nrow = 1)) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    panel.border = element_rect(fill = NA, size = 0.25),
    axis.line = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.title = element_blank(), 
    legend.position = 'bottom',
    legend.background = element_rect(fill = NA),
    legend.margin = margin(c(0,0,0,0)))

gdays + get_legend(gmap) + plot_layout(heights = c(10,1))
ggsave('_figures/s6_sequencedays.png', width = 6, height = 4.5, dpi = 600)

FIG S7: Sequence intensity boxplots

stats_all <- 
  map(.x = 1801:1810, .f = function(huc) {
    rbind(
      df_hist[[which(grid_hucrep[]==huc)]] %>% 
        filter(seq) %>% 
        filter(year(date) %in% 1981:2010) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(source = 'hist'),
      df_ssp245[[which(grid_hucrep[]==huc)]] %>% 
        filter(seq) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(source = 'ssp245'),
      df_ssp585[[which(grid_hucrep[]==huc)]] %>% 
        filter(seq) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(source = 'ssp585')) %>% 
    mutate(
      source_name = case_when(
        source == 'hist' ~ 'Historic',
        source == 'ssp245' ~ 'SSP 2-4.5',
        source == 'ssp585' ~ 'SSP 5-8.5') %>% 
        factor %>% factor(levels = levels(.)[3:1]),
      huc = huc,
      huc4 = paste0('HUC',huc))
  }) %>% reduce(rbind)

## add panel labels
stats_all <- stats_all %>% 
  group_by(huc) %>% 
  mutate(id = 1:length(huc)) %>% 
  ungroup %>% 
  mutate(tag = ifelse(id==1, paste0('(',letters[huc-1800],')'), NA))
ggplot() + 
  # stat_boxplot_custom(
  #   aes(x = maxivt, y = source_name, group = paste(source,ens), fill = source),
  #   qs = c(.05, .25, .5, .75, .95),
  #   outlier.size = 0.5, size = 0.25, outlier.alpha = 1, alpha = 0.5,
  #   show.legend = FALSE) + 
  geom_boxplot(
    data = stats_all %>%
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens,huc4) %>%
      summarize(
        y05 = quantile(maxivt,0.05),
        y25 = quantile(maxivt,0.25),
        y50 = quantile(maxivt,0.50),
        y75 = quantile(maxivt,0.75),
        y95 = quantile(maxivt,0.95)),
    aes(
      xmin = y05, xlower = y25, xmiddle = y50, xupper = y75, xmax = y95,
      y = source_name, group = paste(source_name,ens), fill = source_name),
    stat = 'identity', position = 'dodge',
    width = 0.8, size = 0.25, alpha = 0.5, show.legend = FALSE) +
  geom_point(
    data = stats_all %>% 
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens,huc4) %>% 
      mutate(
        y05 = quantile(maxivt,0.05),
        y95 = quantile(maxivt,0.95)) %>% 
      ungroup %>% group_by(source) %>% 
      mutate(
        dodge = seq(-0.32,0.32,length.out = length(unique(ens)))[ens]) %>% 
      ungroup %>% 
      filter(maxivt < y05 | maxivt > y95),
    aes(x = maxivt, y = as.numeric(factor(source_name))+dodge),
    size = 0.5) + 
  geom_text(
    data = stats_all,
    aes(x = 1650, y = 3.25, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') + 
  facet_wrap(~huc4, ncol = 2, dir = 'v', strip.position = 'left') +
  scale_fill_manual(values = scico(7, palette = 'roma')) +
  scale_x_continuous('Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10)) +
  geom_hline(yintercept = 0.4) + 
  ggtitle('Sequence Intensity') +
  theme(
    text = element_text(family = 'Segoe UI', size = 10),
    panel.grid.major.x = element_line(size = 0.25),
    strip.background = element_rect(color = NA, fill = 'grey95'),
    strip.text = element_text(margin = margin(2,2,2,2)),
    strip.placement = 'inside',
    axis.title.y = element_blank(),
    axis.line = element_line(color = NA)) 
ggsave('_figures/s7_intensity.png', width = 6, height = 7, dpi = 600)

FIG S8: Sequence duration boxplots

ggplot(stats_all) + 
  # stat_boxplot_custom(
  #   aes(x = duration, y = source_name, group = paste(source,ens), fill = source),
  #   qs = c(.05, .25, .5, .75, .95),
  #   outlier.size = 0.5, size = 0.25, outlier.alpha = 1, alpha = 0.5,
  #   show.legend = FALSE) + 
  geom_boxplot(
    data = stats_all %>%
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens,huc4) %>%
      summarize(
        y05 = quantile(duration,0.05),
        y25 = quantile(duration,0.25),
        y50 = quantile(duration,0.50),
        y75 = quantile(duration,0.75),
        y95 = quantile(duration,0.95)),
    aes(
      xmin = y05, xlower = y25, xmiddle = y50, xupper = y75, xmax = y95,
      y = source_name, group = paste(source_name,ens), fill = source_name),
    stat = 'identity', position = 'dodge',
    width = 0.8, size = 0.25, alpha = 0.5, show.legend = FALSE) +
  geom_point(
    data = stats_all %>% 
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens,huc4) %>% 
      mutate(
        y05 = quantile(duration,0.05),
        y95 = quantile(duration,0.95)) %>% 
      ungroup %>% group_by(source) %>% 
      mutate(
        dodge = seq(-0.32,0.32,length.out = length(unique(ens)))[ens]) %>% 
      ungroup %>% 
      filter(duration < y05 | duration > y95),
    aes(x = duration, y = as.numeric(factor(source_name))+dodge),
    size = 0.5) + 
  geom_text(
    aes(x = 140, y = 3.25, label = tag),
    family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
  facet_wrap(~huc4, ncol = 2, dir = 'v', strip.position = 'left') +
  scale_fill_manual(values = scico(7, palette = 'roma')[7:5]) +
  scale_x_continuous('Duration (days)', breaks = 30*(0:10)) +
  geom_hline(yintercept = 0.4) + 
  ggtitle('Sequence Duration') +
  theme(
    text = element_text(family = 'Segoe UI', size = 10),
    panel.grid.major.x = element_line(size = 0.25),
    strip.background = element_rect(color = NA, fill = 'grey95'),
    strip.text = element_text(margin = margin(2,2,2,2)),
    strip.placement = 'inside',
    axis.title.y = element_blank(),
    axis.line = element_line(color = NA)) 
ggsave('_figures/s8_duration.png', width = 6, height = 7, dpi = 600)