library(readr)
library(dplyr)
library(ggplot2)
library(lubridate)
library(sf)
library(mapview)
library(patchwork)
library(here)

source(here("src/analysis/global_analysis_functions.R"))

Overview

This is the velocity analysis for the Drought paper. The weekly velocity records are calculated from instantaneous (15-min) records at four USGS flow stations. Stations include Cache at Ryer Island Ferry, San Joaquin River at Jersey Point, Middle River, and Old River at Bacon Island (see map below). The weekly velocity data is combined with Delta Outflow from Dayflow and drought classifications.

vel_weekly <- readRDS(here("data/processed/hydrology/velocity.rds"))

lat_long <- read_sf(here("data/processed/spatial/velocity_coords.shp"))
# Prepare weekly velocity data for plots
vel_weekly_c <- vel_weekly %>% 
  mutate(
    # Log transform Outflow
    Outflow_log = log(Outflow),
    # Rename stations for more descriptive facet labels
    Station = case_match(
      Station,
      "Cache" ~ "Cache Slough (USGS stations\n11455350 and 11455385)",
      "MDM" ~ "Middle River at Middle River, CA\n(USGS station 11312676)",
      "OLD" ~ "Old River at Bacon Island, CA\n(USGS station 11313405)",
      "SJJ" ~ "SJR at Jersey Point, CA\n(USGS station 11337190)"
    ),
    # Apply custom order to YearType and Station
    YearType = factor(YearType, levels = c("Critical", "Dry", "Below Normal", "Above Normal", "Wet")),
    Station = factor(
      Station,
      levels = c(
        "Cache Slough (USGS stations\n11455350 and 11455385)",
        "SJR at Jersey Point, CA\n(USGS station 11337190)",
        "Middle River at Middle River, CA\n(USGS station 11312676)",
        "Old River at Bacon Island, CA\n(USGS station 11313405)"
      )
    )
  )

Station map

Here is a map of the selected USGS stations flow stations in the Delta - from north to south the stations include Cache Slough (USGS 11455350 - inactive RYI station, 11455385 - active RYF station), Jersey Point (11337190), Old River at Bacon Island (11313405), and Middle River (11312676)

mapview(lat_long, layer.name = "Velocity", size = 10, label = "site")

Figures

Velocity and Delta Outflow by station

Next up is exploring the relationship with delta outflow, and mean net velocity and maximum absolute tidal velocity at four USGS flow stations. Net velocities only increase considerably during periods of high outflow in Cache Slough during wet years. Higher net velocities in Cache Slough reflect periods of Yolo Bypass inundation - when water management activities divert floodwaters from the mainstem Sacramento River. Tidal velocities generally fall within a similar range across all water year types but exhibit a downward trend during periods of high outflow. The negative velocity direction is consistent with periods of higher freshwater flow that mute tidal forcing.

plt_mean_net_vel_vs_outflow <- vel_weekly_c %>%
  ggplot(aes(x = Outflow_log, y = MeanNetVel, color = YearType)) +
  color_pal_yrtype(aes_type = "color", scale_title = "Water Year\nHydrologic\nClassification") +
  geom_point(alpha = 0.4, size = 2) +
  facet_wrap(vars(Station), scales = "free") +
  geom_hline(yintercept = 0, linetype = "solid", col = "black") +
  labs(
    x = expression(log(Outflow)~(m^{3}~s^{-1})),
    y = expression(Mean~Net~Velocity~(m~s^{-1}))
  ) +
  theme_bw()

plt_mean_net_vel_vs_outflow

plt_max_abs_tidal_vel_vs_outflow <- vel_weekly_c %>%
  ggplot(aes(x = Outflow_log, y = MaxAbsTidalVel, color = YearType, shape = TideSign)) +
  geom_point(alpha = 0.4, size = 2) +
  color_pal_yrtype(aes_type = "color", scale_title = "Water Year\nHydrologic\nClassification") +
  facet_wrap(vars(Station), scales = "free") +
  labs(
    x = expression(log(Outflow)~(m^{3}~s^{-1})),
    y = expression(Maximum~Absolute~Tidal~Velocity~(m~s^{-1})),
    shape = "Sign"
  ) +
  theme_bw()

plt_max_abs_tidal_vel_vs_outflow

# Combine figures into one multipanel figure
plt_vel_vs_outflow <- plt_mean_net_vel_vs_outflow / plt_max_abs_tidal_vel_vs_outflow + 
  plot_annotation(tag_levels = "A") +
  plot_layout(guides = "collect")

Delta Outflow

To start, delta outflow is similar at the start of most water years except for wet years when outflow is greater at the beginning. Smoothing the data of each water year type shows outflow ‘peaks’ earlier in the water year during critical and drier years relative to wetter water years. This pattern likely reflects water management activities to increase reservoir storage.

plt_outflow_wyday <- vel_weekly_c %>%
  ggplot(aes(x = WYday, y = Outflow_log, color = YearType)) +
  geom_smooth() +
  color_pal_yrtype(aes_type = "color", scale_title = "Water Year\nHydrologic\nClassification") +
  scale_x_continuous(
    name = "Day of Water Year",
    breaks = seq.int(0, 360, by = 60),
    expand = expansion(mult = 0.005)
  ) +
  ylab(expression(log(Outflow)~(m^{3}~s^{-1}))) +
  theme_bw()

plt_outflow_wyday

Velocity with WY index across stations

Across all water years, mean net velocity is most always positive at Cache and Jersey but only positive in the wet seasons of wet years at the Middle and Old River stations. For much of the year, the export pumps create net landward velocities in Middle and Old Rivers that is demonstrated by negative velocities. Only during the wet season of wet years does increased freshwater flow overcome the effects of the pumps and create net positive velocities.

plt_mean_net_vel_wyday <- vel_weekly_c %>%
  ggplot(aes(x = WYday, y = MeanNetVel, color = YearType)) +
  geom_smooth() +
  color_pal_yrtype(aes_type = "color", scale_title = "Water Year\nHydrologic\nClassification") +
  facet_wrap(vars(Station)) +
  geom_hline(yintercept = 0, linetype = "solid", col = "black") +
  scale_x_continuous(
    name = "Day of Water Year",
    breaks = seq.int(0, 360, by = 60),
    expand = expansion(mult = 0.02)
  ) +
  ylab(expression(Mean~Net~Velocity~(m~s^{-1}))) +
  theme_bw()

plt_mean_net_vel_wyday

plt_max_abs_tidal_vel_wyday <- vel_weekly_c %>%
  ggplot(aes(x = WYday, y = MaxAbsTidalVel, color = YearType)) +
  geom_smooth() +
  color_pal_yrtype(aes_type = "color", scale_title = "Water Year\nHydrologic\nClassification") +
  facet_wrap(vars(Station), scales = "free") +
  scale_x_continuous(
    name = "Day of Water Year",
    breaks = seq.int(0, 360, by = 60),
    expand = expansion(mult = 0.02)
  ) +
  ylab(expression(Maximum~Absolute~Tidal~Velocity~(m~s^{-1}))) +
  theme_bw()

plt_max_abs_tidal_vel_wyday

# Combine figures into one multipanel figure
plt_vel_wyday <- plt_mean_net_vel_wyday / plt_max_abs_tidal_vel_wyday + 
  plot_annotation(tag_levels = "A") +
  plot_layout(guides = "collect")

Export Figures

# Define file path for manuscript figures
fp_plots <- here("results/figures")

# Velocity vs Outflow Figure
ggsave(
  file.path(fp_plots, "vel_vs_outflow.jpg"),
  plot = plt_vel_vs_outflow,
  width = 7,
  height = 9,
  units = "in",
  dpi = 300
)

# Delta Outflow by WY day of year Figure
ggsave(
  file.path(fp_plots, "outflow_wyday.jpg"),
  plot = plt_outflow_wyday,
  width = 5,
  height = 3,
  units = "in",
  dpi = 300
)

# Velocity by WY day of year Figure
ggsave(
  file.path(fp_plots, "vel_wyday.jpg"),
  plot = plt_vel_wyday,
  width = 7,
  height = 9,
  units = "in",
  dpi = 300
)