This markdown script reproduces the figures and numerical results 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), excepting Figure 2. Please note that the figures in the markdown file may not exactly match the published version because of formatting and scale differences.

Setup


Load packages & functions

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

## additional packages not in setup.R
library(patchwork)
library(pammtools)
library(assertthat)
library(devtools)
library(lubridate)

## progress bar
progress <- FALSE
if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)

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

## download custom boxplot function for Figures XX & XX
source_url('https://github.com/BoulderCodeHub/CRSSIO/blob/master/R/stat_boxplot_custom.R?raw=TRUE')

Load MERRA-2 catalogs

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

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

## date 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)

2) Methodology

2.2.1) Sequence Identification Algorithm

FIG 1: Example sequence timeseries


Start-end parameter for example sequence timeseries:

startend <- 
  df_3hr[[167]] %>%
  filter(wateryear(ts) %in% 1981:2010) %>% 
  pull(ivt) %>% median(na.rm = TRUE)

cat(paste('Example start-end parameter value:', round(startend,1), 'kg/m/s'))
## Example start-end parameter value: 91.9 kg/m/s
## pull example data from historic record
df_example <- df_3hr[[167]] %>% 
  filter(wateryear(ts) == 2016) %>%
  filter(month(ts) %in% 10:12) %>% 
  filter(ts >= ymd_hms('2015-10-7 00:00:00')) %>%
  mutate(threshold1 = cbind(rolling5, startend) %>% apply(1, min)) %>% 
  mutate(threshold2 = ifelse(seq.maxrolling > 250, rolling5, threshold1))
  
## create timeseries plot
g <- ggplot(df_example) + 
  geom_hline(yintercept = 250, color = 'grey70') +
  geom_text(
    data = data.frame(x = max(df_example$ts), y = 250),
    aes(x=x, y=y, label = '250 kg/m/s'),
    # nudge_x = 3.1e5, nudge_y = 30, 
    hjust = 1, vjust = -0.5, nudge_x = 8.5e5,
    family = 'Segoe UI', size = 8/.pt,fontface = 'italic', color = 'grey30') + 
  geom_hline(yintercept = startend, color = 'grey70', linetype = 'dashed') +
  geom_text(
    data = data.frame(x = max(df_example$ts), y = 92),
    aes(x=x, y=y, label = '92 kg/m/s'),
    # nudge_x = 3.6e5, nudge_y = 30, 
    hjust = 1, vjust = -0.5, nudge_x = 8.5e5,
    family = 'Segoe UI', size = 8/.pt, fontface = 'italic', color = 'grey30') + 
  geom_ribbon(
    aes(x = ts, ymin = threshold1, ymax = rolling5),
    fill = roma.colors[4], alpha = 0.25) +
  geom_ribbon(
    aes(x = ts, ymin = threshold1, ymax = threshold2),
    fill = roma.colors[4], alpha = 0.75) +
  geom_line(aes(x = ts, y = ivt, color = '3-Hour IVT')) + 
  geom_line(aes(x = ts, y = rolling5, color = '5-Day Moving Average IVT'), size = 0.75) + 
  scale_color_manual('', values = c('grey50', 'black')) + 
  scale_x_datetime(
    'Time (weeks)', date_breaks = '1 week', 
    limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
    expand = expansion(mult = c(0.025, 0.01))) + 
  scale_y_origin('IVT (kg/m/s)', breaks = 250*(0:10)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    legend.text = element_text(size = 8),
    axis.text.x = element_blank(), 
    legend.position = c(0.1,0.8), legend.background = element_blank(),
    panel.grid.major.y = element_line(size = 0.25))

## create bar chart of AR days
g.ar <- ggplot(df_example) + 
  pammtools::geom_stepribbon(
    aes(x = ts, ymin = 0, ymax = as.numeric(ar)), 
    fill = roma.colors[5]) + 
  geom_hline(yintercept = 0, color = roma.colors[5]) +
  geom_text(
    data = data.frame(x = min(df_example$ts), y = 0),
    aes(x=x, y=y, label = 'AR'),
    hjust = 0, vjust = -0.5,
    family = 'Segoe UI', size = 8/.pt, fontface = 'plain') + 
  scale_x_datetime(
    limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
    expand = expansion(mult = c(0.025, 0.01))) + 
  theme_void()

## create bar chart of sequence days
g.seq <- ggplot(df_example) + 
  pammtools::geom_stepribbon(
    aes(x = ts, ymin = 0, ymax = as.numeric(seq)), 
    fill = roma.colors[4]) + 
  geom_hline(yintercept = 0, color = roma.colors[4]) +
  geom_hline(yintercept = 1.5, color = 'white') + 
  geom_text(
    data = data.frame(x = min(df_example$ts), y = 0),
    aes(x=x, y=y, label = 'Sequence'),
    hjust = 0, vjust = -0.5,
    family = 'Segoe UI', size = 8/.pt, fontface = 'plain') + 
  scale_x_datetime(
    limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
    expand = expansion(mult = c(0.025, 0.01))) + 
  theme_void()

## combine into one plot and save
ggarrange(g, g.ar, g.seq, nrow = 3, heights = c(9,1,1.5), align = 'v')
ggsave('_figures/fig1_sequence.png', width = 6, height = 2, units = 'in', dpi = 600)

4) Results: Sequence Characteristics

4.1) Historical Characteristics

Numerical results: Annual frequency statistics

Statewide average annual sequence frequency:

nseq_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_hist[[i]])) NA else {
      temp <- df_hist[[i]] %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        filter(seq) %>% 
        group_by(ens, seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop')
      nrow(temp)/length(unique(temp$ens))
    }
  }
nseq_ref <- nseq_ref / 30
if (progress) cat('\n')

cat(paste(
  'Mean annual frequency:', 
  Mean(nseq_ref) %>% round(2), 'sequences/year \n'))
cat(paste(
  'Standard deviation:', 
  sd(nseq_ref, na.rm = TRUE) %>% round(2), 'sequences/year \n'))
## Mean annual frequency: 1.47 sequences/year 
## Standard deviation: 1.06 sequences/year


Statewide percentage of years with no sequences:

nzero <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_hist[[i]])) NA else {
      temp <- df_hist[[i]] %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        filter(seq) %>% 
        group_by(ens, seq.count) %>% 
        summarize(wy = wateryear(date[1]), duration = seq.duration[1], .groups = 'drop') %>% 
        count(wy,ens) %>% 
        right_join(
          expand.grid(wy = 1981:2010, ens = unique(temp$ens)), 
          by = c('wy','ens')) %>% 
        mutate(n = setNA(n,0))
      temp %>% 
        group_by(ens) %>% 
        summarize(zero = sum(n==0)) %>% 
        pull(zero) %>% mean
    }
  }
if (progress) cat('\n')

cat(paste(
  'Mean number of years with no sequences:', 
  Mean(nzero) %>% round(1), '\n'))
cat(paste(
  'Mean percentage of years with no sequences:', 
  Mean(nzero/30) %>% percent(accuracy = 0.1)))
## Mean number of years with no sequences: 9.7 
## Mean percentage of years with no sequences: 32.5%


Statewide average ARs per sequence in the historic record:

nar_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_hist[[i]])) {
      df_hist[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        group_by(ens,seq.count) %>%
        summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>% 
        summarize(nar.mean = Mean(nar)) %>% unlist %>% unname
    } else NA
  }
if (progress) cat('\n')

cat(paste(
  'Mean number of ARs per sequence:', 
  Mean(nar_ref) %>% round(1), '\n'))
cat(paste(
  'Standard deviation of ARs per sequence:', 
  sd(nar_ref, na.rm = TRUE) %>% round(2), '\n'))
## Mean number of ARs per sequence: 2.4 
## Standard deviation of ARs per sequence: 0.24

FIG 3: Historical Sequence Intensity and Duration

historical <-
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (i %in% index_ca) {
      df_hist[[i]] %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        filter(seq) %>% 
        group_by(seq.count) %>% 
        summarize(
          maxivt = seq.maxivt[1], duration = seq.duration[1], ens = ens[1]) %>% 
        summarize(
          maxivt = mean(maxivt), duration = mean(duration), 
          n.seq = length(seq.count)/length(unique(ens))/30) %>% 
        cbind(
          df_hist[[i]] %>%
            filter(wateryear(date) %in% 1981:2010) %>% 
            filter(ar) %>% 
            group_by(ar.count) %>% 
            summarize(pct.ar = mean(seq)) %>% 
            summarize(pct.ar = mean(pct.ar)))
    } else rep(NA,4)
  }
historical <- historical %>% 
  cbind(raster.df(grid_ca) %>% select(-value)) %>% 
  slice(index_ca)
g1 <- ggplot(historical) + 
  geom_point(aes(x=x, y=y, color = maxivt, size = n.seq)) + 
  geom_sf(data = california, fill = NA, size = 0.25) + 
  geom_point(aes(x=x, y=y, color = maxivt, size = n.seq), alpha = 0.75) + 
  scale_color_gradient(
    'Average \nPeak 24-hr IVT \n(kg/m/s)',
    low = 'white', high = 'darkred') +
  scale_size_area(max_size = 4, guide = guide_none()) + 
  theme(
    legend.position = c(0.92, 0.75),
    axis.title = element_blank(),
    axis.text.x = element_blank(),
    # axis.text.x = element_text(angle = 90, vjust = 0.5),
    text = element_text(family = 'Segoe UI', size = 8))

g2 <- ggplot(historical) + 
  geom_point(aes(x=x, y=y, color = duration, size = n.seq)) + 
  geom_sf(data = california, fill = NA, size = 0.25) + 
  geom_point(aes(x=x, y=y, color = duration, size = n.seq), alpha = 0.75) + 
  scale_color_distiller(
    'Average \nDuration \n(days)',
    palette = 'Blues', direction = 1, breaks = 4*(0:10)) +
  scale_size_area(max_size = 4, guide = guide_none()) + 
  theme(
    legend.position = c(0.85, 0.75),
    axis.title = element_blank(),
    # axis.text.y = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    text = element_text(family = 'Segoe UI', size = 8))

g3 <- ggplot(historical) + 
  geom_point(aes(x=x, y=y, size = n.seq)) + 
  scale_size_area('Annual Sequence Frequency', max_size = 4, breaks = 1:3) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    legend.position = 'bottom')
  
g1 + g2 + plot_spacer() + 
  plot_layout(nrow = 3, heights = c(10,10,1)) +
  inset_element(
    get_legend(g3), 
    left = 0.1, bottom = 0, top = 1, right = 1, 
    ignore_tag = TRUE, align_to = 'full', on_top = TRUE) + 
  plot_annotation(tag_levels = 'a', tag_prefix = '(', tag_suffix = ')') &
  theme(
    plot.tag = element_text(family = 'Segoe UI', face = 'bold', size = 10))
ggsave('_figures/fig3_historic.png', width = 3, height = 6, units = 'in', dpi = 600)

4.2) Projected Frequency Change

Numerical results: Projected change in annual frequency


Create dataframes

nseq_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(seq) %>% 
        group_by(ens, seq.count, decade) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop') %>% 
        count(ens, decade) %>% 
        group_by(decade) %>% 
        summarize(nn = sum(n)/length(unique(ens))) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        mutate(nn = setNA(nn,0)) %>% 
        pull(nn)
    }
  }
nseq_ssp245 <- nseq_ssp245 / 10
if (progress) cat('\n')

nseq_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(seq) %>% 
        group_by(ens, seq.count, decade) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop') %>% 
        count(ens, decade) %>% 
        group_by(decade) %>% 
        summarize(nn = sum(n)/length(unique(ens))) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        mutate(nn = setNA(nn,0)) %>% 
        pull(nn)
    }
  }
nseq_ssp585 <- nseq_ssp585 / 10
if (progress) cat('\n')

nseq_ssp245_eoc <-
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp245[[i]])) NA else {
      temp <- df_ssp245[[i]] %>%
        filter(wateryear(date) %in% 2061:2090) %>% 
        filter(seq) %>%
        group_by(ens, seq.count) %>%
        summarize(duration = seq.duration[1], .groups = 'drop')
      nrow(temp)/length(unique(temp$ens))
    }
  }
nseq_ssp245_eoc <- nseq_ssp245_eoc / 30
if (progress) cat('\n')

nseq_ssp585_eoc <-
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp585[[i]])) NA else {
      temp <- df_ssp585[[i]] %>%
        filter(wateryear(date) %in% 2061:2090) %>% 
        filter(seq) %>%
        group_by(ens, seq.count) %>%
        summarize(duration = seq.duration[1], .groups = 'drop')
      nrow(temp)/length(unique(temp$ens))
    }
  }
nseq_ssp585_eoc <- nseq_ssp585_eoc / 30
if (progress) cat('\n')


Change in sequence frequency by decade:

frequency <- 
  left_join(
    nseq_ssp245 %>% as.data.frame %>%
      cbind(hist = nseq_ref) %>% 
      cbind(huc = grid_hucrep[]) %>% 
      filter(huc > 0) %>% 
      pivot_longer(cols = starts_with('V'), names_to = 'decade') %>% 
      rename(ssp245 = value),
    nseq_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 (Medium)',
      sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>% 
      factor %>% factor(levels = levels(.))) %>%
  mutate(
    id = 1:nrow(.),
    tag = ifelse(id %in% 1:2, paste0('(', letters[id], ')'), NA))

frequency %>% 
  group_by(decade_name, sim) %>% 
  summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>% 
  pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
  gt() %>% 
  tab_header('Change in Sequence Annual Frequency') %>% 
  tab_spanner('SSP2-4.5', columns = ends_with('245')) %>% 
  tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
  fmt_number(starts_with('diff'), decimals = 2) %>% 
  cols_label(
    decade_name = 'Decade',
    diff.mean.ssp245 = 'Mean',
    diff.sd.ssp245 = 'Std.Dev.',
    diff.mean.ssp585 = 'Mean',
    diff.sd.ssp585 = 'Std.Dev') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Change in Sequence Annual Frequency
Decade SSP2-4.5 SSP5-8.5
Mean Std.Dev. Mean Std.Dev
2020s 0.13 0.11 0.25 0.20
2030s 0.37 0.10 0.25 0.19
2040s 0.22 0.19 0.54 0.15
2050s 0.51 0.15 0.82 0.27
2060s 0.59 0.19 0.78 0.14
2070s 0.59 0.08 0.86 0.31
2080s 0.72 0.16 1.13 0.17

FIG 4: Projected change in sequence annual frequency by decade

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

g.additive <- 
  ggplot(frequency) + 
  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.5) + 
  geom_point(
    aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)), 
    size = 0.75) + 
  geom_text(
    aes(x = 1.5, y = 1.6, 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(
    # 'Frequency Change Relative to WY 1981-2010',
    'Change in Sequences/Year',
    # labels = percent, 
    expand = expansion(mult = 0.1)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    panel.grid.major.y = element_line(size = 0.25),
    strip.background = element_rect(color = NA, fill = 'grey95'),
    strip.text = element_text(size = 8, margin = margin(2,0,2,0)),
    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.title = element_text(margin = margin(0,0,0,0)),
    legend.text = element_text(margin = margin(0,0,0,0)),
    legend.margin = margin(c(0,0,0,0)))
gmap_blank <- gmap + guides(fill = guide_none())

layout = 'AAAB\nCCCC'
g.additive + gmap_blank + get_legend(gmap) + plot_spacer() + 
  plot_layout(design = layout, heights = c(10,1))
ggsave('_figures/fig4.png', width = 6, height = 3.5, units = 'in', dpi = 600)

Change in number of ARs per sequence

nar_ssp245 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp245[[i]])) {
      df_ssp245[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(seq) %>% 
        group_by(ens, decade, seq.count) %>% 
        summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>% 
        group_by(decade) %>% 
        summarize(nar = mean(nar)) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        pull(nar)
    } else rep(NA,7)
  }
if (progress) cat('\n')

nar_ssp585 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp585[[i]])) {
      df_ssp585[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(seq) %>% 
        group_by(ens, decade, seq.count) %>% 
        summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>% 
        group_by(decade) %>% 
        summarize(nar = mean(nar)) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        pull(nar)
    } else rep(NA,7)
  }
if (progress) cat('\n')
nar <- 
  left_join(
    nar_ssp245 %>% as.data.frame %>%
      cbind(hist = nar_ref) %>% 
      cbind(huc = grid_hucrep[]) %>% 
      filter(huc > 0) %>% 
      pivot_longer(cols = starts_with('V'), names_to = 'decade') %>% 
      rename(ssp245 = value),
    nar_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 (Medium)',
      sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>% 
      factor %>% factor(levels = levels(.)))

nar %>% 
  group_by(decade_name, sim) %>% 
  summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>% 
  pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
  gt() %>% 
  tab_header('Change in ARs per Sequence') %>% 
  tab_spanner('SSP2-4.5', columns = ends_with('245')) %>% 
  tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
  fmt_number(starts_with('diff'), decimals = 2) %>% 
  cols_label(
    decade_name = 'Decade',
    diff.mean.ssp245 = 'Mean',
    diff.sd.ssp245 = 'Std.Dev.',
    diff.mean.ssp585 = 'Mean',
    diff.sd.ssp585 = 'Std.Dev') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Change in ARs per Sequence
Decade SSP2-4.5 SSP5-8.5
Mean Std.Dev. Mean Std.Dev
2020s −0.02 0.24 0.07 0.20
2030s 0.06 0.24 0.21 0.18
2040s 0.00 0.11 0.25 0.24
2050s 0.08 0.19 0.21 0.21
2060s 0.29 0.21 0.46 0.28
2070s 0.24 0.29 0.44 0.22
2080s 0.32 0.12 0.59 0.27

Change in percentage of ARs falling within sequences

pctar_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_hist[[i]])) {
      df_hist[[i]] %>% 
        filter(ar) %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        group_by(ens,ar.count) %>%
        summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>% 
        pull(pctar) %>% mean
    } else NA
  }
if (progress) cat('\n')

pctar_ssp245 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp245[[i]])) {
      df_ssp245[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(ar) %>% 
        group_by(ens, decade, ar.count) %>% 
        summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>% 
        group_by(decade) %>% 
        summarize(pctar = mean(pctar)) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        pull(pctar)
    } else rep(NA,7)
  }
if (progress) cat('\n')

pctar_ssp585 <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp585[[i]])) {
      df_ssp585[[i]] %>% 
        mutate(yr = year(date)) %>% 
        right_join(decades, by = 'yr') %>% 
        filter(ar) %>% 
        group_by(ens, decade, ar.count) %>% 
        summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>% 
        group_by(decade) %>% 
        summarize(pctar = mean(pctar)) %>% 
        right_join(data.frame(decade = 1:7), by = 'decade') %>% 
        arrange(decade) %>% 
        pull(pctar)
    } else rep(NA,7)
  }
if (progress) cat('\n')
pctar <- 
  left_join(
    pctar_ssp245 %>% as.data.frame %>%
      cbind(hist = pctar_ref) %>% 
      cbind(huc = grid_hucrep[]) %>% 
      filter(huc > 0) %>% 
      pivot_longer(cols = starts_with('V'), names_to = 'decade') %>% 
      rename(ssp245 = value),
    pctar_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 (Medium)',
      sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>% 
      factor %>% factor(levels = levels(.)))

pctar %>% 
  group_by(decade_name, sim) %>% 
  summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>% 
  pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
  gt() %>% 
  tab_header('Change in Percentage of ARs Within Sequences') %>% 
  tab_spanner('SSP2-4.5', columns = ends_with('245')) %>% 
  tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
  fmt_percent(starts_with('diff'), decimals = 1) %>% 
  cols_label(
    decade_name = 'Decade',
    diff.mean.ssp245 = 'Mean',
    diff.sd.ssp245 = 'Std.Dev.',
    diff.mean.ssp585 = 'Mean',
    diff.sd.ssp585 = 'Std.Dev') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Change in Percentage of ARs Within Sequences
Decade SSP2-4.5 SSP5-8.5
Mean Std.Dev. Mean Std.Dev
2020s −4.2% 7.7% −1.3% 3.5%
2030s 1.4% 5.8% −1.2% 5.8%
2040s −1.9% 10.3% 2.5% 8.8%
2050s 1.1% 7.5% 5.7% 5.4%
2060s −0.9% 8.8% 1.5% 5.6%
2070s 1.1% 8.3% 2.2% 8.4%
2080s 3.3% 8.1% 8.2% 4.6%

4.3) Projected Intensity Change

Numerical results: Projected relative change in sequence intensity


Create dataframes

stats_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_hist[[i]])) {
      df_hist[[i]] %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        filter(seq) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        summarize(
          maxivt.med = median(maxivt),
          maxivt.q95 = quantile(maxivt,0.95),
          duration.med = median(duration),
          duration.q95 = quantile(duration,0.95))
    } else rep(NA,4)
  }
if (progress) cat('\n')

stats_ssp245_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp245[[i]])) {
      df_ssp245[[i]] %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        filter(seq) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        summarize(
          maxivt.med = median(maxivt, na.rm = TRUE),
          maxivt.q95 = quantile(maxivt,0.95, na.rm = TRUE),
          duration.med = median(duration, na.rm = TRUE),
          duration.q95 = quantile(duration,0.95, na.rm = TRUE))
    } else rep(NA,4)
  }
if (progress) cat('\n')

stats_ssp585_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp585[[i]])) {
      df_ssp585[[i]] %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        filter(seq) %>% 
        group_by(ens,seq.count) %>%
        summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
        summarize(
          maxivt.med = median(maxivt, na.rm = TRUE),
          maxivt.q95 = quantile(maxivt,0.95, na.rm = TRUE),
          duration.med = median(duration, na.rm = TRUE),
          duration.q95 = quantile(duration,0.95, na.rm = TRUE))
    } else rep(NA,4)
  }
if (progress) cat('\n')

stats <- 
  rbind(
    df_hist[[which(grid_hucrep[]==1804)]] %>% 
      filter(wateryear(date) %in% 1981:2010) %>% 
      filter(seq) %>% 
      group_by(ens,seq.count) %>%
      summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>% 
      mutate(source = 'hist'),
    df_ssp245[[which(grid_hucrep[]==1804)]] %>% 
      filter(wateryear(date) %in% 2061:2090) %>% 
      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[]==1804)]] %>% 
      filter(wateryear(date) %in% 2061:2090) %>% 
      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]))

Filter distribution statistics to focus on sequence intensity

intensity <- 
  cbind(
    (stats_ssp245_eoc / stats_ref) %>%
      select(starts_with('maxivt')) %>%
      setNames(gsub('maxivt', 'ssp245', names(.))),
    (stats_ssp585_eoc / stats_ref) %>%
      select(starts_with('maxivt')) %>%
      setNames(gsub('maxivt', 'ssp585', names(.)))) %>%
    cbind(raster.df(grid_ca) %>% select(-value)) %>%
    slice(index_ca) %>%
    pivot_longer(cols = c(-x,-y)) %>%
    separate(name, into = c('sim', 'stat')) %>%
    mutate(
      sim_name = case_when(
        sim == 'ssp245' ~ 'SSP 2-4.5',
        sim == 'ssp585' ~ 'SSP 5-8.5') %>% 
        factor %>% factor(levels = levels(.)[2:1]),
      stat_name = case_when(
        stat == 'med' ~ 'Median',
        stat == 'q95' ~ '95th Percentile') %>% 
        factor %>% factor(levels = levels(.)[3:1]))


Relative intensity change for HUC1804:

left_join(
  stats_ssp245_eoc %>% 
    select(starts_with('maxivt')) %>% 
    cbind(huc = grid_hucrep[]) %>% 
    filter(huc == 1804) %>% 
    select(-huc) %>% 
    pivot_longer(everything(), names_to = 'stat', values_to = 'ssp245'),
  stats_ref %>% 
    select(starts_with('maxivt')) %>% 
    cbind(huc = grid_hucrep[]) %>% 
    filter(huc == 1804) %>% 
    select(-huc) %>% 
    pivot_longer(everything(), names_to = 'stat', values_to = 'hist'), 
  by = 'stat') %>% 
  mutate(
    stat = gsub('maxivt.', '', stat),
    stat_name = case_when(
      stat == 'med' ~ 'Median',
      stat == 'q95' ~ '95th Percentile') %>% 
      factor %>% factor(levels = levels(.)[3:1])) %>% 
  mutate(fact = ssp245/hist, pct = fact-1) %>% 
  select(stat_name, hist, ssp245, fact, pct) %>% 
  gt %>% 
  tab_header(
    'Relative Change Calculation', 
    subtitle = 'Sequence Intensity (Peak 24-Hour IVT, kg/m/s)') %>% 
  fmt_number(c(hist,ssp245), decimals = 0) %>% 
  fmt_number(fact, decimals = 2) %>% 
  fmt_percent(pct, decimals = 0) %>% 
  cols_label(
    stat_name = 'Statistic', 
    hist = 'Historic',
    ssp245 = 'SSP2-4.5', 
    fact = 'Change Factor',
    pct = 'Percent Change') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Relative Change Calculation
Sequence Intensity (Peak 24-Hour IVT, kg/m/s)
Statistic Historic SSP2-4.5 Change Factor Percent Change
Median 539 576 1.07 7%
95th Percentile 715 832 1.16 16%


Statewide average relative intensity change:

intensity %>% 
  group_by(sim, stat) %>% 
  summarize(mean = Mean(value)-1, sd = sd(value, na.rm = TRUE), .groups = 'drop') %>% 
  # mutate(pct.mean = round(100*(fact.mean-1),1), pct.sd = round(100*fact.sd,1)) %>%
  mutate(sim_name = ifelse(sim == 'ssp245', 'SSP2-4.5', 'SSP5-8.5')) %>% 
  select(-sim) %>% 
  pivot_wider(names_from = stat, values_from = c(mean,sd), names_sep = '.') %>% 
  select(sim_name, ends_with('med'), ends_with('iqr'), ends_with('q95')) %>%
  gt %>% 
  tab_header(
    'Relative Change in Statewide Average Statistics',
    subtitle = 'Sequence Intensity (Peak 24-Hour IVT, kg/m/s)') %>% 
  tab_spanner('Median', columns = ends_with('med')) %>% 
  tab_spanner('IQR', columns = ends_with('iqr')) %>% 
  tab_spanner('95th Percentile', columns = ends_with('q95')) %>% 
  fmt_percent(c(starts_with('mean'), starts_with('sd')), decimals = 1) %>% 
  cols_label(
    sim_name = 'Emission Scenario',
    mean.med = 'Mean', sd.med = 'Std.Dev.',
    mean.q95 = 'Mean', sd.q95 = 'Std.Dev.') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Relative Change in Statewide Average Statistics
Sequence Intensity (Peak 24-Hour IVT, kg/m/s)
Emission Scenario Median 95th Percentile
Mean Std.Dev. Mean Std.Dev.
SSP2-4.5 6.5% 4.3% 12.7% 6.2%
SSP5-8.5 11.6% 5.6% 19.1% 5.7%

FIG 5: Sequence intensity change

Top row (a,b):

a <- ggplot(stats) +
  geom_ribbon(
    aes(x = maxivt, group = source, ymin = source_name, ymax = source_name, fill = source),
    color = NA, alpha = 0) +
  # 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.25,
  #   show.legend = FALSE) + 
  geom_boxplot(
    data = stats %>%
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens) %>%
      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),
    stat = 'identity', position = 'dodge', 
    width = 0.8, size = 0.25, alpha = 0.25, show.legend = FALSE) +
  geom_point(
    data = stats %>% 
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens) %>% 
      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) + 
  scale_fill_manual(
    'Climate Scenario', 
    labels = c('Historic', 'SSP 2-4.5', 'SSP 5-8.5'), 
    values = scico(7, palette = 'roma')[3:1], 
    guide = guide_legend(
      override.aes = list(
        color = scico(7, palette = 'roma')[3:1], size = 0.25, alpha = 0.5))) +
  scale_x_continuous(
    'Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10), 
    limits = c(250,1250), expand = expansion(mult = c(0.005,0.05))) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    panel.grid.major.x = element_line(size = 0.25),
    axis.title.y = element_blank(),
    legend.position = 'top',
    legend.margin = margin(1,1,1,1)) + #c(0.9, 0.75)) + 
  ggtitle('Sequence Intensity') + 
  geom_segment(
    x = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
    xend = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
    y = 2.6, yend = 3.4,
    color = scico(11, palette = 'roma')[4]) +
  geom_segment(
    x = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
    xend = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
    y = 1.6, yend = 2.4,
    color = scico(11, palette = 'roma')[3]) + 
  geom_segment(
    x = stats %>% filter(source == 'ssp585') %>% pull(maxivt) %>% median,
    xend = stats %>% filter(source == 'ssp585') %>% pull(maxivt) %>% median,
    y = 0.6, yend = 1.4,
    color = scico(11, palette = 'roma')[1])
a

b <- ggplot() +
  geom_tile(
    data = raster.df(grid_hucrep),
    aes(x=x, y=y, fill = value==1804)) + 
  scale_fill_manual(
    values = c('white', 'grey50'), guide = guide_none()) +
  geom_tile(
    data = raster.df(grid_hucrep*grid_ca/grid_ca) %>% arrange(value==1804),
    aes(x=x, y=y, color = value==1804), fill = NA) + 
  scale_color_manual(
    values = c('grey70', 'black'), na.value = NA, guide = guide_none()) + 
  geom_sf(data = cal, fill = NA, color = 'black', size = 0.25) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.title = element_blank(), 
    legend.background = element_rect(fill = NA))
b

Bottom row (c,d,e):

c <- ggplot(stats %>% filter(source != 'ssp585')) + 
  ## medians
  geom_vline(
    xintercept = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
    size = 0.75, color = baker[3]) +
  geom_vline(
    xintercept = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
    size = 0.75, color = baker[3]) +
  ## boxplot
  # stat_boxplot_custom(
  #   aes(x = maxivt, y = source_name, fill = source),
  #   qs = c(.05, .25, .5, .75, .95),
  #   outlier.size = 0.5, size = 0.25, outlier.alpha = 0, alpha = 0.25,
  #   width = 0.5,
  #   show.legend = FALSE) + 
  geom_boxplot(
    data = stats %>% filter(source != 'ssp585') %>%
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source) %>%
      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), fill = source),
    stat = 'identity', width = 0.5, size = 0.2, alpha = 0.25, show.legend = FALSE) +
  scale_fill_manual(values = scico(7, palette = 'roma')[c(3:1)]) +
  ## 95th quantiles
  geom_vline(
    xintercept = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% quantile(0.95),
    size = 0.75, color = roma.colors[5], linetype = 'dashed') +
  geom_vline(
    xintercept = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% quantile(0.95),
    size = 0.75, color = roma.colors[5], linetype = 'dashed') +
  ## plot specs
  scale_x_continuous(
    'Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10), 
    limits = c(250,1250), expand = expansion(mult = c(0.005,0.05))) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    panel.grid.major.x = element_line(size = 0.25),
    axis.title.x = element_text(hjust = 0.25),
    axis.title.y = element_blank()) 
c

de <- intensity %>% 
  mutate(
    id = 1:nrow(.),
    tag = ifelse(id %in% 1:3, paste0('(', letters[4:6], ')'), NA),
    sim_name = factor(sim_name, levels = rev(levels(sim_name)))) %>% 
  select(-id) %>% 
  filter(stat != 'iqr') %>% 
  ggplot() +
  geom_vline(xintercept = 0, color = 'grey50', size = 0.35) +
  geom_density(
    aes(x = value-1, group = sim_name, fill = sim_name, color = sim_name),
    alpha = 0.35, size = 0.35) +
  scale_color_manual(
    'Climate Scenario', values = scico(7, palette = 'roma')[2:1]) + 
  scale_fill_manual(
    'Climate Scenario', values = scico(7, palette = 'roma')[2:1]) + 
  guides(color = guide_none(), fill = guide_none()) + 
  facet_wrap(~stat_name, nrow = 1) + 
  scale_x_continuous(
    'Percent Change Relative to 1981-2010', labels = percent) +
  scale_y_origin('Probability Density') +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    legend.margin = margin(0,0,0,0),
    legend.position = 'top', #c(0.85,0.5),
    panel.grid.major.x = element_line(size = 0.25),
    strip.placement = 'outside',
    strip.background = element_rect(fill = 'grey95', color = NA),
    strip.text = element_text(size = 8, margin = margin(2,2,2,2)))
de

ggsave(a, filename = '_figures/fig5a.png', width = 4.25, height = 2.25, dpi = 600)
ggsave(b, filename = '_figures/fig5b.png', width = 1.75, height = 1.75, dpi = 600)
ggsave(c, filename = '_figures/fig5c.png', width = 4.25, height = 1.5, dpi = 600)
ggsave(de, filename = '_figures/fig5de.png', width = 3, height = 1.75, dpi = 600)

The individual figure panels were annotated and combined in Inkscape.

4.4) Duration

Numerical results: Projected relative change in sequence duration


Filter distribution statistics to focus on sequence duration

duration <- 
  cbind(
    (stats_ssp245_eoc / stats_ref) %>%
      select(starts_with('duration')) %>%
      setNames(gsub('duration', 'ssp245', names(.))),
    (stats_ssp585_eoc / stats_ref) %>%
      select(starts_with('duration')) %>%
      setNames(gsub('duration', 'ssp585', names(.)))) %>%
    cbind(raster.df(grid_ca) %>% select(-value)) %>%
    slice(index_ca) %>%
    pivot_longer(cols = c(-x,-y)) %>%
    separate(name, into = c('sim', 'stat')) %>%
    mutate(
      sim_name = case_when(
        sim == 'ssp245' ~ 'SSP 2-4.5',
        sim == 'ssp585' ~ 'SSP 5-8.5') %>% 
        factor %>% factor(levels = levels(.)[2:1]),
      stat_name = case_when(
        stat == 'med' ~ 'Median',
        stat == 'iqr' ~ 'Interquartile Range (IQR)', 
        stat == 'q95' ~ '95th Percentile') %>% 
        factor %>% factor(levels = levels(.)[3:1]))


Relative duration change for HUC1804:

left_join(
  stats_ssp245_eoc %>% 
    select(starts_with('duration')) %>% 
    cbind(huc = grid_hucrep[]) %>% 
    filter(huc == 1804) %>% 
    select(-huc) %>% 
    pivot_longer(everything(), names_to = 'stat', values_to = 'ssp245'),
  stats_ref %>% 
    select(starts_with('duration')) %>% 
    cbind(huc = grid_hucrep[]) %>% 
    filter(huc == 1804) %>% 
    select(-huc) %>% 
    pivot_longer(everything(), names_to = 'stat', values_to = 'hist'), 
  by = 'stat') %>% 
  mutate(
    stat = gsub('duration.', '', stat),
    stat_name = case_when(
      stat == 'med' ~ 'Median',
      stat == 'q95' ~ '95th Percentile') %>% 
      factor %>% factor(levels = levels(.)[3:1])) %>% 
  mutate(fact = ssp245/hist, pct = fact-1) %>% 
  select(stat_name, hist, ssp245, fact, pct) %>% 
  gt %>% 
  tab_header(
    'Relative Change Calculation', 
    subtitle = 'Sequence Duration (days)') %>% 
  fmt_number(c(hist,ssp245), decimals = 0) %>% 
  fmt_number(fact, decimals = 2) %>% 
  fmt_percent(pct, decimals = 0) %>% 
  cols_label(
    stat_name = 'Statistic', 
    hist = 'Historic',
    ssp245 = 'SSP2-4.5', 
    fact = 'Change Factor',
    pct = 'Percent Change') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Relative Change Calculation
Sequence Duration (days)
Statistic Historic SSP2-4.5 Change Factor Percent Change
Median 20 23 1.15 15%
95th Percentile 52 53 1.02 2%


Statewide average relative duration change:

duration %>% 
  group_by(sim, stat) %>% 
  summarize(mean = Mean(value)-1, sd = sd(value, na.rm = TRUE), .groups = 'drop') %>% 
  mutate(sim_name = ifelse(sim == 'ssp245', 'SSP2-4.5', 'SSP5-8.5')) %>% 
  pivot_wider(names_from = stat, values_from = c(mean,sd), names_sep = '.') %>% 
  select(sim_name, ends_with('med'), ends_with('iqr'), ends_with('q95')) %>%
  gt %>% 
  tab_header(
    'Relative Change in Statewide Average Statistics',
    subtitle = 'Sequence Duration (days)') %>% 
  tab_spanner('Median', columns = ends_with('med')) %>% 
  tab_spanner('IQR', columns = ends_with('iqr')) %>% 
  tab_spanner('95th Percentile', columns = ends_with('q95')) %>% 
  fmt_percent(c(starts_with('mean'), starts_with('sd')), decimals = 1) %>% 
  cols_label(
    sim_name = 'Emission Scenario',
    mean.med = 'Mean', sd.med = 'Std.Dev.',
    mean.q95 = 'Mean', sd.q95 = 'Std.Dev.') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Relative Change in Statewide Average Statistics
Sequence Duration (days)
Emission Scenario Median 95th Percentile
Mean Std.Dev. Mean Std.Dev.
SSP2-4.5 10.9% 10.5% 7.6% 12.4%
SSP5-8.5 17.9% 10.0% 28.4% 16.8%

FIG 6: Sequence duration change

a <- ggplot(stats) + 
  geom_ribbon(
    aes(x = duration, group = source, ymin = source_name, ymax = source_name, fill = source),
    color = NA, alpha = 0) + 
  # 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.35,
  #   show.legend = FALSE) + 
  geom_boxplot(
    data = stats %>%
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens) %>%
      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),
    stat = 'identity', position = 'dodge',
    width = 0.8, size = 0.25, alpha = 0.35, show.legend = FALSE) +
  geom_point(
    data = stats %>% 
      mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>% 
      group_by(source_name,source,ens) %>% 
      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) + 
  scale_fill_manual(
    'Climate Scenario', 
    labels = c('Historic', 'SSP 2-4.5', 'SSP 5-8.5'), 
    values = scico(7, palette = 'roma')[5:7], 
    guide = guide_legend(
      override.aes = list(
        color = scico(7, palette = 'roma')[5:7], size = 0.25, alpha = 0.5))) +
  scale_x_continuous('Duration (days)', breaks = 30*(0:10)) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    panel.grid.major.x = element_line(size = 0.25),
    axis.title.y = element_blank(),
    legend.position = 'top',
    legend.margin = margin(1,1,1,1)) + 
  ggtitle('Sequence Duration') + 
  geom_segment(
    x = stats %>% filter(source == 'hist') %>% pull(duration) %>% median,
    xend = stats %>% filter(source == 'hist') %>% pull(duration) %>% median,
    y = 2.6, yend = 3.4,
    color = scico(11, palette = 'roma')[8]) +
  geom_segment(
    x = stats %>% filter(source == 'ssp245') %>% pull(duration) %>% median,
    xend = stats %>% filter(source == 'ssp245') %>% pull(duration) %>% median,
    y = 1.6, yend = 2.4,
    color = scico(11, palette = 'roma')[9]) + 
  geom_segment(
    x = stats %>% filter(source == 'ssp585') %>% pull(duration) %>% median,
    xend = stats %>% filter(source == 'ssp585') %>% pull(duration) %>% median,
    y = 0.6, yend = 1.4,
    color = scico(11, palette = 'roma')[11])
a

bc <- duration %>% 
  mutate(
    id = 1:nrow(.),
    tag = ifelse(id %in% 1:3, paste0('(', letters[2:4], ')'), NA),
    sim_name = factor(sim_name, levels = rev(levels(sim_name)))) %>% 
  select(-id) %>% 
  filter(stat != 'iqr') %>% 
  ggplot() +
  geom_vline(xintercept = 0, color = 'grey50', size = 0.35) +
  geom_density(
    aes(x = value-1, group = sim_name, fill = sim_name, color = sim_name),
    alpha = 0.35, size = 0.35) +
  scale_color_manual('Climate\nScenario', values = scico(7, palette = 'roma')[6:7]) + 
  scale_fill_manual('Climate\nScenario', values = scico(7, palette = 'roma')[6:7]) + 
  guides(color = guide_none(), fill = guide_none()) + 
  facet_wrap(~stat_name, nrow = 1) + 
  scale_x_continuous(
    'Percent Change Relative to 1981-2010', labels = percent) +
  scale_y_origin('Probability Density\n') +
  coord_cartesian(xlim = c(NA,1.25)) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    legend.margin = margin(0,0,0,0),
    legend.position = c(0.9,0.5),
    panel.grid.major.x = element_line(size = 0.25),
    strip.placement = 'outside',
    strip.background = element_rect(fill = 'grey95', color = NA),
    strip.text = element_text(margin = margin(2,2,2,2)), 
    plot.margin = margin(5,5,5,5))
bc

ggsave(a, filename = '_figures/fig6a.png', width = 4.25, height = 2.25, dpi = 600)
ggsave(bc, filename = '_figures/fig6bc.png', width = 4.18, height = 1.75, dpi = 600)

The individual figure panels were annotated and combined in Inkscape.

Numerical results: Super-sequences


Create dataframes

extremes_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_hist[[i]])) rep(NA,5) else {
      temp <- df_hist[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop')
      temp %>% 
        group_by(ens) %>% 
        summarize(n = sum(duration>60)/30) %>% 
        right_join(data.frame(ens = unique(df_hist[[i]]$ens))) %>%
        mutate(n = setNA(n,0)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% 
        pull(n)
    }
  } %>% apply(1, Max)
if (progress) cat('\n')

extremes_ssp245_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp245[[i]])) rep(NA,5) else {
      temp <- df_ssp245[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop')
      temp %>% 
        group_by(ens) %>% 
        summarize(n = sum(duration>60)/30) %>% 
        right_join(data.frame(ens = unique(df_ssp245[[i]]$ens))) %>%
        mutate(n = setNA(n,0)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% 
        pull(n)
    }
  } %>% apply(1, Max)
if (progress) cat('\n')

extremes_ssp585_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (is.null(df_ssp585[[i]])) rep(NA,5) else {
      temp <- df_ssp585[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop')
      temp %>% 
        group_by(ens) %>% 
        summarize(n = sum(duration>60)/30) %>% 
        right_join(data.frame(ens = unique(df_ssp585[[i]]$ens))) %>%
        mutate(n = setNA(n,0)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% 
        pull(n)
    }
  } %>% apply(1, Max)
if (progress) cat('\n')


Historical vs. projected future spatial coverage of super-sequences:

rbind(
  c('Historic', 'single', Sum(extremes_ref>0) / sum(!is.na(extremes_ref))),
  c('SSP2-4.5', 'single', Sum(extremes_ssp245_eoc>0) / sum(!is.na(extremes_ssp245_eoc))),
  c('SSP5-8.5', 'single', Sum(extremes_ssp585_eoc>0) / sum(!is.na(extremes_ssp585_eoc))),
  c('Historic', 'multiple', Sum(extremes_ref>(1/30)) / sum(!is.na(extremes_ref))),
  c('SSP2-4.5', 'multiple', Sum(extremes_ssp245_eoc>(1/30)) / sum(!is.na(extremes_ssp245_eoc))),
  c('SSP5-8.5', 'multiple', Sum(extremes_ssp585_eoc>(1/30)) / sum(!is.na(extremes_ssp585_eoc)))) %>% 
  as.data.frame %>% setNames(c('source', 'super', 'pct')) %>% 
  mutate(pct = toNumber(pct)) %>% 
  pivot_wider(names_from = super, values_from = pct) %>% 
  gt %>% 
  tab_header(
    'Super-Sequence Coverage', 
    subtitle = 'Percent of Grid Cells Affected over 30 Years (1981-2010 or 2021-2050)') %>% 
  fmt_percent(c(single, multiple), decimals = 1) %>% 
  cols_label(
    source = 'Emission Scenario',
    single = 'Cells with >0 Super-Sequences',
    multiple = 'Cells with Multiple Super-Sequences') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Super-Sequence Coverage
Percent of Grid Cells Affected over 30 Years (1981-2010 or 2021-2050)
Emission Scenario Cells with >0 Super-Sequences Cells with Multiple Super-Sequences
Historic 73.3% 39.4%
SSP2-4.5 86.7% 67.3%
SSP5-8.5 98.8% 93.9%

FIG 7: Sequence annual frequency as a function of duration

bucket_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_hist[[i]])) {
      temp <- df_hist[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 1981:2010) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
      (summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>% 
        as.matrix %>% t %>% as.data.frame
    } else rep(NA,6)
  }
if (progress) cat('\n')
bucket_ssp245_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp245[[i]])) {
      temp <- df_ssp245[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
      (summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>% 
        as.matrix %>% t %>% as.data.frame
    } else rep(NA,6)
  }
if (progress) cat('\n')
bucket_ssp585_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp585[[i]])) {
      temp <- df_ssp585[[i]] %>% 
        filter(seq) %>% 
        filter(wateryear(date) %in% 2061:2090) %>% 
        group_by(ens,seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop') %>% 
        mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
      (summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>% 
        as.matrix %>% t %>% as.data.frame
    } else rep(NA,6)
  }
if (progress) cat('\n')

bucket <- 
  rbind(
    bucket_ref %>% 
      cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>% 
      slice(index_ca) %>% 
      mutate(id = index_ca, source = 'hist'),
    bucket_ssp245_eoc %>% 
      cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>% 
      slice(index_ca) %>% 
      mutate(id = index_ca, source = 'ssp245'),
    bucket_ssp585_eoc %>% 
      cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>% 
      slice(index_ca) %>% 
      mutate(id = index_ca, source = 'ssp585')) %>% 
  pivot_longer(cols = c(-id,-source,-huc4), names_to = 'bucket', values_to = 'freq') %>% 
  mutate(
    bucket_label = case_when(
      bucket == '(0,14]' ~ 'Duration \u2264 14 Days',
      bucket == '(14,30]' ~ 'Duration = 15-30 Days',
      bucket == '(30,60]' ~ 'Duration = 31-60 Days',
      bucket == '(60,200]' ~ 'Duration > 60 Days') %>% 
      factor(levels = c(
        'Duration \u2264 14 Days', 'Duration = 15-30 Days',
        'Duration = 31-60 Days', 'Duration > 60 Days'))) %>% 
  pivot_wider(names_from = source, values_from = freq)
temp <- bucket %>% 
  filter(bucket_label == levels(bucket_label)[4]) %>% 
  select(hist, ssp245, ssp585) %>% 
  apply(2, Max)
ratio <- (temp/temp[1])[-1]

uppers <- bucket %>% 
  group_by(bucket) %>% 
  summarize(across(c('hist', 'ssp245', 'ssp585'), Max)) %>% 
  select(-bucket)

g.high <- list()
for (i in 1:4) {
  temp <- 
    ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) + 
    geom_point(aes(x = hist, y = ssp245), size = 0.75)
  dx <- ggplot_build(temp)$layout$panel_params[[1]]$x$breaks %>% diff %>% .[1]
  g.high[[i]] <- 
    ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) + 
    geom_point(aes(x = hist, y = ssp245), size = 0.75) + 
    facet_wrap(~bucket_label) + 
    scale_x_origin('Historic (1921-2010) Annual Frequency') + 
    scale_y_origin(
      'SSP2-4.5 Annual\nFrequency (WY 2061-2090)', 
      breaks = seq(0, uppers$ssp245[i]+2*dx, 2*dx)) + 
    geom_parity() + 
    coord_cartesian(
      clip = 'off', xlim = c(0, uppers$hist[i]), ylim = c(0, ratio[1]*1.1*uppers$hist[i])) + 
    theme(
      text = element_text(family = 'Segoe UI', size = 8),
      panel.grid.major = element_line(size = 0.25),
      axis.title = element_blank(),
      strip.background = element_rect(color = NA, fill = 'grey95'),
      strip.text = element_text(margin = margin(2,2,2,2)))
}
g.high[[1]] <- g.high[[1]] + theme(axis.title.y = element_text())

g.low <- list()
for (i in 1:4) {
  temp <- 
    ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) + 
    geom_point(aes(x = hist, y = ssp585), size = 0.75)
  dx <- ggplot_build(temp)$layout$panel_params[[1]]$x$breaks %>% diff %>% .[1]
  g.low[[i]] <- 
    ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) + 
    geom_point(aes(x = hist, y = ssp585), size = 0.75) + 
    facet_wrap(~bucket_label) + 
    scale_x_origin('Historic (1921-2010)\nAnnual Frequency') + 
    scale_y_origin(
      'SSP5-8.5 Annual\nFrequency (WY 2061-2090)', 
      breaks = seq(0, uppers$ssp585[i]+2*dx, 2*dx)) + 
    geom_parity() + 
    coord_cartesian(
      clip = 'off', xlim = c(0, uppers$hist[i]), ylim = c(0, ratio[2]*uppers$hist[i])) + 
    theme(
      text = element_text(family = 'Segoe UI', size = 8),
      panel.grid.major = element_line(size = 0.25),
      axis.title = element_blank(),
      strip.background = element_rect(color = NA, fill = 'grey95'),
      strip.text = element_text(margin = margin(2,2,2,2)))
}
g.low[[1]] <- g.low[[1]] + theme(axis.title.y = element_text())

layout <- "ABCD\nEFGH"
g.high[[1]] + g.high[[2]] + g.high[[3]] + g.high[[4]] + #g.high[[5]] + 
  g.low[[1]] + g.low[[2]] + g.low[[3]] + g.low[[4]] + #g.low[[5]] +  
  plot_layout(
    design = layout, 
    widths = rep(1,5), heights = c(3,5)) + 
  plot_annotation(
    caption = 'Historic Annual Frequency (WY 1981-2010)',
    tag_levels = list(letters[c(1,3,5,7,2,4,6,8)]), tag_prefix = '(', tag_suffix = ')', ) &
  theme(
    plot.tag = element_text(
      family = 'Segoe UI', face = 'bold', size = 10, margin = margin(-5,-5,-5,-5)),
    plot.caption = element_text(size = 8, hjust = 0.5, margin = margin(0,0,0,0)))
ggsave('_figures/fig7_returnperiod.png', width = 6, height = 3.5, units = 'in', dpi = 600)

5) Discussion


Create dataframes

exceed_ref <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_hist[[i]])) {
      df_hist[[i]] %>% 
        mutate(wy = wateryear(date)) %>% 
        filter(wy %in% 1981:2010) %>% 
        filter(seq) %>% 
        group_by(ens, wy, seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
        summarize(exceed60 = (sum(duration>60)/length(duration)) > 0, .) %>% 
        right_join(expand.grid(ens = unique(df_hist[[i]]$ens), wy = 1981:2010), by = c('ens', 'wy')) %>%
        mutate(exceed60 = setNA(exceed60, FALSE)) %>% 
        summarize(across(starts_with('exceed'), mean)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% pull(exceed60)
    } else rep(NA,5)
  }
if (progress) cat('\n')

exceed_ssp245_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp245[[i]])) {
      df_ssp245[[i]] %>% 
        mutate(wy = wateryear(date)) %>% 
        filter(wy %in% 2061:2090) %>% 
        filter(seq) %>% 
        group_by(ens, wy, seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
        summarize(exceed60 = (sum(duration>60)/length(duration)) > 0) %>% 
        right_join(expand.grid(ens = unique(df_ssp245[[i]]$ens), wy = 2061:2090), by = c('ens', 'wy')) %>%
        mutate(exceed60 = setNA(exceed60, FALSE)) %>% 
        summarize(across(starts_with('exceed'), mean)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% pull(exceed60)
    } else rep(NA,4)
  }
if (progress) cat('\n')

exceed_ssp585_eoc <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (!is.null(df_ssp585[[i]])) {
      df_ssp585[[i]] %>% 
        mutate(wy = wateryear(date)) %>% 
        filter(wy %in% 2061:2090) %>% 
        filter(seq) %>% 
        group_by(ens, wy, seq.count) %>% 
        summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
        summarize(exceed60 = (sum(duration>60)/length(duration)) > 0) %>% 
        right_join(expand.grid(ens = unique(df_ssp585[[i]]$ens), wy = 2061:2090), by = c('ens', 'wy')) %>%
        mutate(exceed60 = setNA(exceed60, FALSE)) %>% 
        summarize(across(starts_with('exceed'), mean)) %>% 
        right_join(data.frame(ens = 1:5), by = 'ens') %>% 
        arrange(ens) %>% pull(exceed60)
    } else rep(NA,4)
  }
if (progress) cat('\n')

Percentage of years with super-sequences:

data.frame(
  `Historic` = exceed_ref %>% apply(1, Mean),
  `SSP2-4.5` = exceed_ssp245_eoc %>% apply(1, Mean),
  `SSP5-8.5` = exceed_ssp585_eoc %>% apply(1, Mean)) %>%
  slice(index_ca) %>% 
  pivot_longer(everything()) %>%
  group_by(name) %>%
  # as_tibble() %>% 
  summarize(across(value, .fns = list(mean = mean, sd = sd, rp = function(x) 1/mean(x)))) %>% 
  # pivot_longer(-name, names_to = 'stat')
  gt %>% 
  tab_header('Percent of Years with Super-Sequences') %>% 
  fmt_markdown(name) %>% 
  fmt_percent(c(value_mean, value_sd), decimals = 1) %>% 
  fmt_number(value_rp, decimals = 1) %>% 
  cols_label(
    name = 'Emission Scenario',
    value_mean = 'Mean (%)',
    value_sd = 'Std. Dev. (%)',
    value_rp = 'Return Period (yrs)') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')
Percent of Years with Super-Sequences
Emission Scenario Mean (%) Std. Dev. (%) Return Period (yrs)

Historic

1.9% 1.9% 53.0

SSP2.4.5

4.2% 3.6% 23.9

SSP5.8.5

10.5% 6.3% 9.5