## ============================================================
## PACKAGES
## ============================================================
# install.packages("dataRetrieval") # run once if needed
# install.packages("trend") # run once if needed
library(dataRetrieval)
library(trend)
library(knitr)
library(kableExtra)
## ============================================================
## CONFIGURATION (EDIT THIS BLOCK FOR EACH NEW FIRE)
## ============================================================
## Use when:
## - switch to a new fire
## - use different USGS gages
## - plug in a different precip CSV
## ------------------------------------------------------------
## --- USGS site numbers ---
# Fire-affected stream (downstream of burn)
site_stream <- "07103800" # Eagle Creek (fire-affected)
# Control stream (not burned, similar region)
site_control <- "07104000" # Rio Ruidoso (control)
## --- Basin Drainage Areas
basin_stream_area_km2 <- 38
basin_control_area_km2 <-526
# Only change the km2 areas section, the m2 will auto update below
stream_area_m2 <- basin_stream_area_km2 * 1e6
control_area_m2 <- basin_control_area_km2 * 1e6
## --- Names for plots and labels (for titles & legends) ---
fire_name <- "WC Fire (2012)"
basin_stream_name <- "West Monument Creek"
basin_control_name <- "Monument Creek"
## --- Analysis period (calendar dates) ---
start_date <- "2005-01-01"
end_date <- "2018-12-31"
## --- Fire water year split ---
## WY < fire_wy as "pre-fire"
## WY >= fire_wy as "post-fire"
fire_wy <- 2012
## --- Water-year analysis window .
wy_min <- 2006
wy_max <- 2018
## --- Precipitation data (Google Sheets CSV) ---
## This is a regional daily precipitation time series (one station), which we:
precip_url <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vTUUfrM0KaKiFCMEIsD3e68Vsxo_GEno0dX0MbM7ahD0uUwvXJ87XIFbAR_hMySViHxS0T3lXiaM2Ik/pub?output=csv"
## ============================================================
## CONVERSIONS – GLOBAL UNIT CONSTANTS
## ============================================================
## Purpose:
## Centralize all unit conversions used anywhere in the script.
## ------------------------------------------------------------
## --- Discharge conversions ---
cfs_to_m3s <- 0.0283168466 # cubic feet per second → m³/s
seconds_per_day <- 86400
cfs_to_m3day <- cfs_to_m3s * seconds_per_day # 1 cfs → m³/day (≈ 2446.576)
## --- Precipitation conversions ---
inch_to_m <- 0.0254 # inches → meters
inch_to_mm <- 25.4 # inches → millimeters
## --- Area conversions ---
km2_to_m2 <- 1e6 # km² → m²
## --- flow depth conversions ---
m3_to_mm <- 1000 # convert m³ per m² → mm depth
mm_to_m <- 1e-3 # mm → m (rarely needed)
## Now any section in the script can use:
## discharge_m3day <- Q_cfs * cfs_to_m3day
## precip_m <- precip_in * inch_to_m
## area_m2 <- area_km2 * km2_to_m2
## etc.
## ============================================================
## FUNCTIONS – GENERIC & REUSABLE
## ============================================================
## ------------------------------------------------------------
## water_year(date_vec)
## ------------------------------------------------------------
water_year <- function(date_vec) {
yr <- as.integer(format(date_vec, "%Y"))
mo <- as.integer(format(date_vec, "%m"))
wy <- ifelse(mo >= 10, yr + 1, yr)
return(wy)
}
## ------------------------------------------------------------
## add_wy_and_dQ(df)
##
## Input:
## df: raw daily discharge from readNWISdv()
##
## Output:
## df with added columns:
## - date : Date class
## - Q_cfs : discharge (cfs)
## - wy : water year
## - dQ : absolute daily change in discharge
##
## Significance:
## dQ is the day-to-day change in flow. Summed over a water
## year and normalized by total flow, this becomes the
## Richards–Baker Index (RBI), a metric of "flashiness":
## * high RBI → more spiky, flashy hydrograph
## * low RBI → smoother, more baseflow-dominated behavior
## ------------------------------------------------------------
add_wy_and_dQ <- function(df) {
# Convert to Date object and numeric discharge
df$date <- as.Date(df$Date)
df$Q_cfs <- as.numeric(df$X_00060_00003)
# Sort rows by date just in case
df <- df[order(df$date), ]
# Add water year
df$wy <- water_year(df$date)
# Lagged discharge (yesterday’s flow)
Q_lag <- c(NA, head(df$Q_cfs, -1))
# Daily change in discharge (absolute value)
df$dQ <- abs(df$Q_cfs - Q_lag)
return(df)
}
## ------------------------------------------------------------
## compute_flashiness(df, min_days = 300)
##
## Input:
## df: data frame with columns wy, Q_cfs, dQ
##
## Output:
## data.frame with one row per water year:
## - wy : water year
## - flash : Richards–Baker Flashiness Index
## - meanQ : mean daily discharge that WY
## - n_days : how many valid Q days in that WY
##
## Formula:
## flashiness = sum(|dQ|) / sum(Q)
##
## Significance:
## This compresses the entire annual hydrograph into a single
## metric that tells you how "spiky" the flow regime is. Post-
## fire
## * higher flashiness (more event-driven flows)
## * changes in mean Q (e.g., more runoff vs infiltration)
##
## Drop water years with < min_days of data to avoid
## interpreting incomplete records as a "real" change.
## ------------------------------------------------------------
compute_flashiness <- function(df, min_days = 300) {
wy_list <- sort(unique(df$wy))
out <- data.frame(
wy = wy_list,
flash = NA_real_,
meanQ = NA_real_,
n_days = NA_integer_
)
for (i in seq_along(wy_list)) {
w <- wy_list[i]
sub <- df[df$wy == w, ]
n <- sum(!is.na(sub$Q_cfs))
out$n_days[i] <- n
if (n >= min_days) {
out$flash[i] <- sum(sub$dQ, na.rm = TRUE) / sum(sub$Q_cfs, na.rm = TRUE)
out$meanQ[i] <- mean(sub$Q_cfs, na.rm = TRUE)
}
}
# Remove water years that didn't meet min_days requirement
out <- out[!is.na(out$flash), ]
return(out)
}
# --- Helper: run Pettitt on a series ordered by WY and report the WY ---
run_pettitt <- function(wy_vec, value_vec, label = "series") {
ord <- order(wy_vec)
wy_ord <- wy_vec[ord]
x_ord <- value_vec[ord]
pt <- pettitt.test(x_ord)
# change-point index
k <- pt$estimate
cp_wy <- wy_ord[k]
cat("\nPettitt test for", label, "\n")
cat(" Change-point WY:", cp_wy, "\n")
cat(" K statistic :", pt$statistic, "\n")
cat(" p-value :", pt$p.value, "\n\n")
invisible(list(result = pt, cp_wy = cp_wy))
}
## ------------------------------------------------------------
## assign_prepost_colors(period_vec, pre_col, post_col)
##
## Simple helper to consistently map:
## period == "pre" → pre_col
## period == "post" → post_col
## ------------------------------------------------------------
assign_prepost_colors <- function(period_vec, pre_col = "blue", post_col = "orange") {
ifelse(period_vec == "pre", pre_col, post_col)
}
## ------------------------------------------------------------
## Helper for Wilcoxon Rank-Sum Tests (Mann-Whitney)
##
## Inputs:
## df : data frame with at least {response, period}
## response : string, name of column to test (e.g., "flash")
## label : string printed in output header
##
## Output:
## A list containing the test object and the median diff
##
## Prints:
## - variable tested
## - group medians
## - Hodges–Lehmann median difference
## - 95% confidence interval
## - p-value
## ------------------------------------------------------------
## Small helper to build a summary row-set for any MWW test
mww_stats <- function(df, response) {
y <- df[[response]]
grp <- df$period
med_pre <- median(y[grp == "pre"], na.rm = TRUE)
med_post <- median(y[grp == "post"], na.rm = TRUE)
wt <- wilcox.test(
y ~ grp,
exact = FALSE,
conf.int = TRUE,
conf.level = 0.95
)
data.frame(
Statistic = c(
"Pre-fire median",
"Post-fire median",
"Median difference (post - pre)",
"95% CI lower",
"95% CI upper",
"p-value"
),
Value = c(
med_pre,
med_post,
unname(wt$estimate),
wt$conf.int[1],
wt$conf.int[2],
wt$p.value
)
)
}
run_mww_pretty <- function(df, response, label = "Wilcoxon Test") {
# Extract variables
y <- df[[response]]
grp <- df$period
# Compute group medians
med_pre <- median(y[grp == "pre"], na.rm = TRUE)
med_post <- median(y[grp == "post"], na.rm = TRUE)
# Run test with CI + effect size
wt <- wilcox.test(
y ~ grp,
exact = FALSE,
conf.int = TRUE,
conf.level = 0.95
)
# Pretty printed output
cat("\n------------------------------------------------------------\n")
cat(label, "\n")
cat("------------------------------------------------------------\n")
cat("Pre-fire median: ", round(med_pre, 4), "\n")
cat("Post-fire median: ", round(med_post, 4), "\n")
cat("Median difference (post - pre): ", round(wt$estimate, 4), "\n")
cat("95% CI: [", round(wt$conf.int[1], 4), ", ",
round(wt$conf.int[2], 4), "]\n", sep = "")
cat("p-value: ", format.pval(wt$p.value), "\n")
invisible(wt)
}
## ============================================================
## DATA RETRIEVAL & PROCESSING – DISCHARGE
## ============================================================
## Steps:
## 1) Download daily discharge for Eagle Creek (fire) and
## Ruidoso (control) using USGS web services.
## 2) Add water year and daily |dQ|.
## 3) Aggregate by water year to get RBI and mean Q.
## 4) Restrict to water-year window of interest.
##
## Hydrologic significance:
## * Comparing fire-affected vs control basins helps separate
## climate-driven changes from fire-driven changes.
## ------------------------------------------------------------
## 1) Download daily discharge for FIRE (stream) and CONTROL
# Fire-affected site (Eagle Creek)
Q_stream_raw <- readNWISdv(
siteNumbers = site_stream,
parameterCd = "00060", # discharge
startDate = start_date,
endDate = end_date
)
## Warning in readNWISdv(siteNumbers = site_stream, parameterCd = "00060", : NWIS
## servers are slated for decommission. Please begin to migrate to
## read_waterdata_daily.
## GET: https://waterservices.usgs.gov/nwis/dv/?site=07103800&format=waterml%2C1.1&ParameterCd=00060&StatCd=00003&startDT=2005-01-01&endDT=2018-12-31
# Control site (Ruidoso)
Q_control_raw <- readNWISdv(
siteNumbers = site_control,
parameterCd = "00060",
startDate = start_date,
endDate = end_date
)
## Warning in readNWISdv(siteNumbers = site_control, parameterCd = "00060", : NWIS
## servers are slated for decommission. Please begin to migrate to
## read_waterdata_daily.
## GET: https://waterservices.usgs.gov/nwis/dv/?site=07104000&format=waterml%2C1.1&ParameterCd=00060&StatCd=00003&startDT=2005-01-01&endDT=2018-12-31
# Add water year and daily dQ for both
Q_stream <- add_wy_and_dQ(Q_stream_raw)
Q_control <- add_wy_and_dQ(Q_control_raw)
# Compute annual flashiness for each (using all available dates within start/end)
flash_stream_all <- compute_flashiness(Q_stream)
flash_control_all <- compute_flashiness(Q_control)
# Restrict flashiness to chosen WY window (full water years only)
flash_stream <- flash_stream_all[
flash_stream_all$wy >= wy_min & flash_stream_all$wy <= wy_max,
]
flash_control <- flash_control_all[
flash_control_all$wy >= wy_min & flash_control_all$wy <= wy_max,
]
# Also create WY-limited daily discharge for plotting (nice to match precip)
Q_stream_plot <- Q_stream[Q_stream$wy >= wy_min & Q_stream$wy <= wy_max, ]
Q_control_plot <- Q_control[Q_control$wy >= wy_min & Q_control$wy <= wy_max, ]
## ============================================================
## DATA RETRIEVAL & PROCESSING – PRECIPITATION
## ============================================================
## Goal:
## Convert a daily precip station time series into:
## (1) Daily values aligned with WYs
## (2) Water-year totals (precip_wy)
##
## Significance:
## Need this to:
## * Check if post-fire changes in flow could simply be due
## to wetter/drier years.
## * Build Q/P and RBI/P ratios that normalize out climate
## forcing as best as possible.
## ------------------------------------------------------------
## 2) Read and process precipitation from Google Sheets URL
# ===== Load CSV directly from Google Sheets =====
precip_raw <- read.csv(precip_url, stringsAsFactors = FALSE)
# ===== Standardize column names (assumes 2 columns) =====
names(precip_raw) <- c("date_raw", "precip_raw")
# ===== Convert date to proper Date class =====
precip_raw$date <- as.Date(precip_raw$date_raw, format = "%m/%d/%y")
# ===== Clean precipitation values =====
precip_raw$precip_in <- as.character(precip_raw$precip_raw)
# M = missing → NA
precip_raw$precip_in[precip_raw$precip_in == "M"] <- NA
# T = trace → 0
precip_raw$precip_in[precip_raw$precip_in == "T"] <- "0"
# Convert to numeric
precip_raw$precip_in <- as.numeric(precip_raw$precip_in)
# ===== Add water year (Oct–Sep) =====
precip_raw$water_year <- water_year(precip_raw$date)
# ===== Final cleaned dataframe (keep only what you need) =====
precip_df <- precip_raw[, c("date", "precip_in", "water_year")]
names(precip_df)[names(precip_df) == "water_year"] <- "wy"
# Restrict precip to the same WY window (full water years only)
precip_df <- precip_df[precip_df$wy >= wy_min & precip_df$wy <= wy_max, ]
### precip_df now contains:
# date (Date)
# precip_in (numeric, with T=0 and M=NA)
# wy (numeric water year), from wy_min to wy_max
# ============================================================
# ANNUAL (WATER-YEAR) PRECIP TOTALS
# ============================================================
## Significance:
## precip_wy is your climate forcing series. If flashiness or
## Q/P changes while precip_wy stays roughly similar, that’s
## evidence for internal watershed change (e.g., wildfire
## impacts) rather than just "it rained more."
# ------------------------------------------------------------
precip_wy <- aggregate(
precip_in ~ wy,
data = precip_df,
FUN = sum,
na.rm = TRUE
)
## ============================================================
## RUNOFF RATIO (RR) – ANNUAL WATER BALANCE METRIC
## ============================================================
## Purpose:
## Compute the annual runoff ratio for each basin:
##
## RR = (annual runoff depth) / (annual precipitation depth)
##
## where:
## - annual runoff depth = (sum of Q over WY) / basin area
## - precipitation depth = WY total precipitation (converted to mm)
##
## Hydrologic meaning:
## - RR expresses how efficiently precipitation becomes streamflow.
## - RR close to 1 = most rainfall exits the basin as runoff
## - RR << 1 = large losses to ET, infiltration, storage
##
## This provides a long-term water-balance perspective that
## complements:
## * RBI (hydrograph shape)
## * Q/P_peak (event-scale sensitivity)
## ============================================================
## --- 1) Compute annual discharge volume for each basin (m³ / WY)
Q_stream_wy_vol <- aggregate(
Q_cfs ~ wy,
data = Q_stream[!is.na(Q_stream$Q_cfs), ],
FUN = function(x) sum(x * cfs_to_m3day)
)
names(Q_stream_wy_vol)[2] <- "Q_vol_stream_m3"
Q_control_wy_vol <- aggregate(
Q_cfs ~ wy,
data = Q_control[!is.na(Q_control$Q_cfs), ],
FUN = function(x) sum(x * cfs_to_m3day)
)
names(Q_control_wy_vol)[2] <- "Q_vol_control_m3"
## --- 2) Convert discharge volume to runoff depth (mm)
Q_stream_wy_vol$runoff_mm <- (Q_stream_wy_vol$Q_vol_stream_m3 / stream_area_m2) * 1000
Q_control_wy_vol$runoff_mm <- (Q_control_wy_vol$Q_vol_control_m3 / control_area_m2) * 1000
## --- 3) Prepare annual precipitation depth (mm)
precip_wy$P_mm <- precip_wy$precip_in * inch_to_mm
## --- 4) Merge Q_depth and P_depth for RR calculation
RR_stream <- merge(Q_stream_wy_vol[, c("wy", "runoff_mm")],
precip_wy[, c("wy", "P_mm")],
by = "wy")
RR_control <- merge(Q_control_wy_vol[, c("wy", "runoff_mm")],
precip_wy[, c("wy", "P_mm")],
by = "wy")
## --- 5) Compute runoff ratio
RR_stream$RR <- RR_stream$runoff_mm / RR_stream$P_mm
RR_control$RR <- RR_control$runoff_mm / RR_control$P_mm
## --- 6) Restrict to complete water years + WY window
# Use the WYs already accepted as complete in the flashiness section
valid_wy_stream <- flash_stream$wy
valid_wy_control <- flash_control$wy
RR_stream <- RR_stream[RR_stream$wy %in% valid_wy_stream , ]
RR_control <- RR_control[RR_control$wy %in% valid_wy_control , ]
## ============================================================
## Data proceessing: PEAK-FLOW SENSITIVITY (Q_peak / P_peak)
## ============================================================
## Goal:
## Relate the annual maximum daily discharge (Q_peak) to the
## annual maximum daily precipitation (P_peak) in each WY.
##
## Metric:
## Q_peak / P_peak (units: cfs per inch)
## ------------------------------------------------------------
## ---- Annual peak discharge (per WY) ----
Q_stream_peak <- aggregate(
Q_cfs ~ wy,
data = Q_stream[!is.na(Q_stream$Q_cfs), ],
FUN = max
)
names(Q_stream_peak)[2] <- "Q_peak_stream"
Q_control_peak <- aggregate(
Q_cfs ~ wy,
data = Q_control[!is.na(Q_control$Q_cfs), ],
FUN = max
)
names(Q_control_peak)[2] <- "Q_peak_control"
## ---- Annual peak daily precip (per WY) ----
precip_peak <- aggregate(
precip_in ~ wy,
data = precip_df[!is.na(precip_df$precip_in), ],
FUN = max
)
names(precip_peak)[2] <- "P_peak"
## ---- Merge peaks by WY ----
peak_stream <- merge(Q_stream_peak, precip_peak, by = "wy", all = FALSE)
peak_control <- merge(Q_control_peak, precip_peak, by = "wy", all = FALSE)
## ---- Restrict to analysis WY window ----
peak_stream <- peak_stream[peak_stream$wy >= wy_min & peak_stream$wy <= wy_max, ]
peak_control <- peak_control[peak_control$wy >= wy_min & peak_control$wy <= wy_max, ]
## ---- Compute Q/P ratios (cfs per inch) ----
peak_stream$ratio_QP <- peak_stream$Q_peak_stream / peak_stream$P_peak
peak_control$ratio_QP <- peak_control$Q_peak_control / peak_control$P_peak
## ---- Tag pre vs post fire ----
peak_stream$period <- ifelse(peak_stream$wy < fire_wy, "pre", "post")
peak_control$period <- ifelse(peak_control$wy < fire_wy, "pre", "post")
## ---- Drop any WY with missing ratio ----
peak_stream <- peak_stream[!is.na(peak_stream$ratio_QP), ]
peak_control <- peak_control[!is.na(peak_control$ratio_QP), ]
## ============================================================
## PRE / POST FIRE TAGGING (WATER YEAR–BASED)
## ============================================================
## We split all WY-based series into:
## period = "pre" for WY < fire_wy
## period = "post" for WY >= fire_wy
##
## This is the basic grouping used for:
## * Mann–Whitney tests
## * plotting color schemes
## ------------------------------------------------------------
# Tag pre vs post for flashiness at stream (fire site)
flash_stream$period <- ifelse(flash_stream$wy < fire_wy, "pre", "post")
# Same split for control site to keep comparison aligned
flash_control$period <- ifelse(flash_control$wy < fire_wy, "pre", "post")
# Tag pre vs post in precip WY totals
precip_wy$period <- ifelse(precip_wy$wy < fire_wy, "pre", "post")
## --- 7) Tag pre/post fire for RR
RR_stream$period <- ifelse(RR_stream$wy < fire_wy, "pre", "post")
RR_control$period <- ifelse(RR_control$wy < fire_wy, "pre", "post")
## ============================================================
## PLOTTING – BASIC DIAGNOSTIC PLOTS
## ============================================================
## Purpose:
## Visual inspection of:
## * How flashiness evolves through time
## * Whether daily Q patterns look different post-fire
## * How precip behaves through the same window
## ------------------------------------------------------------
## Pre/post color vectors (using helper)
cols_stream <- assign_prepost_colors(flash_stream$period)
cols_control <- assign_prepost_colors(flash_control$period)
cols_precip <- assign_prepost_colors(precip_wy$period)
## For vertical line in time-series plots:
fire_start_date <- as.Date(paste0(fire_wy - 1, "-10-01")) # 1 Oct of fire WY
## 1) Flashiness – stream (fire basin)
plot(
flash_stream$wy, flash_stream$flash,
type = "b", pch = 16, col = cols_stream,
xlab = "Water Year",
ylab = "Flashiness (RBI)",
main = paste("Flashiness –", basin_stream_name,
"(", site_stream, ")", sep = " ")
)
abline(v = fire_wy, lty = 2, col = "red")
legend(
"topleft",
legend = c("Pre-fire", "Post-fire", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## 2) Flashiness – control basin
plot(
flash_control$wy, flash_control$flash,
type = "b", pch = 16, col = cols_control,
xlab = "Water Year",
ylab = "Flashiness (RBI)",
main = paste("Flashiness –", basin_control_name,
"(", site_control, ")", sep = " ")
)
abline(v = fire_wy, lty = 2, col = "red")
legend(
"topleft",
legend = c("Pre-fire", "Post-fire", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## 3) Daily discharge – stream (limited to WY window)
plot(
Q_stream_plot$date, Q_stream_plot$Q_cfs,
type = "l",
col = "steelblue",
xlab = "Date",
ylab = "Discharge (cfs)",
main = paste("Daily Discharge –", basin_stream_name,
"(", site_stream, ")", sep = " ")
)
abline(v = fire_start_date, col = "red", lty = 2)
legend(
"topleft",
legend = c("Discharge", fire_name),
col = c("steelblue", "red"),
lty = c(1, 2),
bty = "n"
)
## 4) Daily discharge – control (limited to WY window)
plot(
Q_control_plot$date, Q_control_plot$Q_cfs,
type = "l",
col = "darkgreen",
xlab = "Date",
ylab = "Discharge (cfs)",
main = paste("Daily Discharge –", basin_control_name,
"(", site_control, ")", sep = " ")
)
abline(v = fire_start_date, col = "red", lty = 2)
legend(
"topleft",
legend = c("Discharge", fire_name),
col = c("darkgreen", "red"),
lty = c(1, 2),
bty = "n"
)
## 5) Daily precipitation time series (WY window)
plot(
precip_df$date, precip_df$precip_in,
type = "h",
col = "gray40",
xlab = "Date",
ylab = "Precipitation (in/day)",
main = "Daily Precipitation (Regional Station)"
)
abline(v = fire_start_date, col = "red", lty = 2)
legend(
"topleft",
legend = c("Daily precip", paste("Start of WY", fire_wy)),
col = c("gray40", "red"),
lty = c(1, 2),
bty = "n"
)
## 6) Annual (water-year) precipitation totals
plot(
precip_wy$wy, precip_wy$precip_in,
type = "b", pch = 16, col = cols_precip,
xlab = "Water Year",
ylab = "Total Precipitation (in / WY)",
main = "Annual (Water-Year) Precipitation"
)
abline(v = fire_wy, col = "red", lty = 2)
legend(
"topleft",
legend = c("Pre-fire WY", "Post-fire WY", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## --- 7) RR (runoff Ratio)
plot(
RR_stream$wy, RR_stream$RR,
type = "b", pch = 16,
col = ifelse(RR_stream$period == "pre", "blue", "orange"),
xlab = "Water Year",
ylab = "Runoff Ratio (RR)",
main = paste("Runoff Ratio –", basin_stream_name)
)
abline(v = fire_wy, col = "red", lty = 2)
legend(
"topleft",
legend = c("Pre-fire WY", "Post-fire WY", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
plot(
RR_control$wy, RR_control$RR,
type = "b", pch = 16,
col = ifelse(RR_control$period == "pre", "blue", "orange"),
xlab = "Water Year",
ylab = "Runoff Ratio (RR)",
main = paste("Runoff Ratio –", basin_control_name)
)
abline(v = fire_wy, col = "red", lty = 2)
legend(
"bottomright",
legend = c("Pre-fire WY", "Post-fire WY", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## ============================================================
## Plot: Peak-flow sensitivity (Q_peak / P_peak)
## ============================================================
## ---- FIRE site ----
cols_QP_stream <- ifelse(peak_stream$period == "pre", "blue", "orange")
plot(
peak_stream$wy, peak_stream$ratio_QP,
type = "b", pch = 16,
col = cols_QP_stream,
xlab = "Water Year",
ylab = "Q_peak / P_peak (cfs per inch)",
main = paste("Peak Flow Sensitivity –", basin_stream_name)
)
abline(v = fire_wy, col = "red", lty = 2)
legend(
"topleft",
legend = c("Pre-fire WY", "Post-fire WY", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## ---- CONTROL site ----
cols_QP_control <- ifelse(peak_control$period == "pre", "blue", "orange")
plot(
peak_control$wy, peak_control$ratio_QP,
type = "b", pch = 16,
col = cols_QP_control,
xlab = "Water Year",
ylab = "Q_peak / P_peak (cfs per inch)",
main = paste("Peak Flow Sensitivity –", basin_control_name)
)
abline(v = fire_wy, col = "red", lty = 2)
legend(
"topleft",
legend = c("Pre-fire WY", "Post-fire WY", fire_name),
col = c("blue", "orange", "red"),
pch = c(16, 16, NA),
lty = c(NA, NA, 2),
bty = "n"
)
## ============================================================
## RBI vs. PRECIPITATION – BASE CASE
## ============================================================
## Here we relate annual flashiness (RBI) to annual precip.
##
## Questions:
## - Do wetter years tend to be flashier?
## - Does the slope or intercept change post-fire?
##
## We first fit regression lines to:
## * all pre-fire years
## * all post-fire years
## ------------------------------------------------------------
# Merge WY flashiness with WY precipitation
flash_precip <- merge(
flash_stream_all,
precip_wy[, c("wy", "precip_in")],
by = "wy"
)
# Tag pre vs post fire
flash_precip$period <- ifelse(flash_precip$wy < fire_wy, "pre", "post")
# Colors
colvec <- ifelse(flash_precip$period == "pre", "blue", "orange")
# RBI vs Precipitation Plot
plot(
flash_precip$precip_in,
flash_precip$flash,
pch = 19,
col = colvec,
xlab = "Annual Precipitation (inches)",
ylab = "Flashiness (RBI)",
main = "RBI vs Precipitation (Pre vs Post Fire)"
)
# Regression lines
pre_mod <- lm(flash ~ precip_in, data = flash_precip[flash_precip$period == "pre", ])
post_mod <- lm(flash ~ precip_in, data = flash_precip[flash_precip$period == "post", ])
#comment out to remove trendline
#abline(pre_mod, col = "blue", lwd = 2)
#abline(post_mod, col = "orange", lwd = 2, lty = 2)
legend(
"topright",
legend = c("Pre-fire", "Post-fire"),
col = c("blue", "orange"),
pch = 19,
lty = c(1, 2),
bty = "n"
)
## ============================================================
## Statistcal Analysis Section
## ============================================================
## ============================================================
## 1) MANN–WHITNEY (WILCOXON) TESTS – PRE vs POST
## ============================================================
## Compare distributions of flashiness, mean Q, and RR
## before vs after fire.
## ------------------------------------------------------------
# Flashiness – stream (fire site)
w_flash_stream <- run_mww_pretty(
flash_stream, "flash",
label = "Flashiness – FIRE Site (Stream)"
)
##
## ------------------------------------------------------------
## Flashiness – FIRE Site (Stream)
## ------------------------------------------------------------
## Pre-fire median: 0.0844
## Post-fire median: 0.11
## Median difference (post - pre): 0.0253
## 95% CI: [-0.0535, 0.0918]
## p-value: 0.61708
# Flashiness – control site
w_flash_control <- run_mww_pretty(
flash_control, "flash",
label = "Flashiness – CONTROL Site"
)
##
## ------------------------------------------------------------
## Flashiness – CONTROL Site
## ------------------------------------------------------------
## Pre-fire median: 0.2516
## Post-fire median: 0.3132
## Median difference (post - pre): 0.0606
## 95% CI: [-0.0119, 0.1626]
## p-value: 0.074146
# Summary table for flashiness (stream vs control)
flash_fire_stats <- mww_stats(flash_stream, "flash")
flash_ctrl_stats <- mww_stats(flash_control, "flash")
mww_flash_table <- data.frame(
Statistic = flash_fire_stats$Statistic,
FIRE_Flash = flash_fire_stats$Value,
CONTROL_Flash = flash_ctrl_stats$Value
)
knitr::kable(
mww_flash_table,
digits = 4,
caption = "Mann–Whitney tests for flashiness (RBI): FIRE vs CONTROL"
)
Mann–Whitney tests for flashiness (RBI): FIRE vs
CONTROL
| Pre-fire median |
0.0844 |
0.2516 |
| Post-fire median |
0.1100 |
0.3132 |
| Median difference (post - pre) |
0.0253 |
0.0606 |
| 95% CI lower |
-0.0535 |
-0.0119 |
| 95% CI upper |
0.0918 |
0.1626 |
| p-value |
0.6171 |
0.0741 |
# Mean Q – stream (fire site)
w_meanQ_stream <- run_mww_pretty(
flash_stream, "meanQ",
label = "Mean Q – FIRE Site (Stream)"
)
##
## ------------------------------------------------------------
## Mean Q – FIRE Site (Stream)
## ------------------------------------------------------------
## Pre-fire median: 0.9628
## Post-fire median: 4.4491
## Median difference (post - pre): 3.4971
## 95% CI: [1.5812, 5.5092]
## p-value: 0.018416
# Mean Q – control site
w_meanQ_control <- run_mww_pretty(
flash_control, "meanQ",
label = "Mean Q – CONTROL Site"
)
##
## ------------------------------------------------------------
## Mean Q – CONTROL Site
## ------------------------------------------------------------
## Pre-fire median: 32.7533
## Post-fire median: 29.2479
## Median difference (post - pre): 2.7666
## 95% CI: [-11.7614, 36.2593]
## p-value: 0.72098
# Summary table for mean Q (stream vs control)
meanQ_fire_stats <- mww_stats(flash_stream, "meanQ")
meanQ_ctrl_stats <- mww_stats(flash_control, "meanQ")
mww_meanQ_table <- data.frame(
Statistic = meanQ_fire_stats$Statistic,
FIRE_MeanQ = meanQ_fire_stats$Value,
CONTROL_MeanQ = meanQ_ctrl_stats$Value
)
knitr::kable(
mww_meanQ_table,
digits = 4,
caption = "Mann–Whitney tests for mean daily discharge: FIRE vs CONTROL"
)
Mann–Whitney tests for mean daily discharge: FIRE vs
CONTROL
| Pre-fire median |
0.9628 |
32.7533 |
| Post-fire median |
4.4491 |
29.2479 |
| Median difference (post - pre) |
3.4971 |
2.7666 |
| 95% CI lower |
1.5812 |
-11.7614 |
| 95% CI upper |
5.5092 |
36.2593 |
| p-value |
0.0184 |
0.7210 |
# Runoff Ratio – stream (fire site)
w_RR_stream <- run_mww_pretty(
RR_stream, "RR",
label = "Runoff Ratio (RR) – FIRE Site (Stream)"
)
##
## ------------------------------------------------------------
## Runoff Ratio (RR) – FIRE Site (Stream)
## ------------------------------------------------------------
## Pre-fire median: 0.0638
## Post-fire median: 0.2377
## Median difference (post - pre): 0.1765
## 95% CI: [0.0936, 0.2365]
## p-value: 0.008221
# Runoff Ratio – control site
w_RR_control <- run_mww_pretty(
RR_control, "RR",
label = "Runoff Ratio (RR) – CONTROL Site"
)
##
## ------------------------------------------------------------
## Runoff Ratio (RR) – CONTROL Site
## ------------------------------------------------------------
## Pre-fire median: 0.1402
## Post-fire median: 0.1401
## Median difference (post - pre): 0.0045
## 95% CI: [-0.065, 0.0649]
## p-value: 0.83032
# Summary table for RR (stream vs control)
RR_fire_stats <- mww_stats(RR_stream, "RR")
RR_ctrl_stats <- mww_stats(RR_control, "RR")
mww_RR_table <- data.frame(
Statistic = RR_fire_stats$Statistic,
FIRE_RR = RR_fire_stats$Value,
CONTROL_RR = RR_ctrl_stats$Value
)
knitr::kable(
mww_RR_table,
digits = 4,
caption = "Mann–Whitney tests for runoff ratio (RR): FIRE vs CONTROL"
)
Mann–Whitney tests for runoff ratio (RR): FIRE vs
CONTROL
| Pre-fire median |
0.0638 |
0.1402 |
| Post-fire median |
0.2377 |
0.1401 |
| Median difference (post - pre) |
0.1765 |
0.0045 |
| 95% CI lower |
0.0936 |
-0.0650 |
| 95% CI upper |
0.2365 |
0.0649 |
| p-value |
0.0082 |
0.8303 |
## ============================================================
## 2) PEAK-FLOW SENSITIVITY (Q_peak / P_peak)
## ============================================================
## Goal:
## Relate the annual maximum daily discharge (Q_peak) to the
## annual maximum daily precipitation (P_peak) in each WY.
##
## Metric:
## Q_peak / P_peak (units: cfs per inch)
## ------------------------------------------------------------
## ---- Wilcoxon tests on Q/P ratio (pretty output) ----
w_ratio_stream <- run_mww_pretty(
peak_stream,
response = "ratio_QP",
label = "Q_peak / P_peak – FIRE Site (Stream)"
)
##
## ------------------------------------------------------------
## Q_peak / P_peak – FIRE Site (Stream)
## ------------------------------------------------------------
## Pre-fire median: 5.1589
## Post-fire median: 12.5217
## Median difference (post - pre): 8.7056
## 95% CI: [1.9247, 25.0331]
## p-value: 0.018416
w_ratio_control <- run_mww_pretty(
peak_control,
response = "ratio_QP",
label = "Q_peak / P_peak – CONTROL Site"
)
##
## ------------------------------------------------------------
## Q_peak / P_peak – CONTROL Site
## ------------------------------------------------------------
## Pre-fire median: 260.0089
## Post-fire median: 294.0789
## Median difference (post - pre): 36.1176
## 95% CI: [-122.0362, 268.6274]
## p-value: 0.61708
## ---- Summary table for Q_peak / P_peak (stream vs control) ----
QP_fire_stats <- mww_stats(peak_stream, "ratio_QP")
QP_ctrl_stats <- mww_stats(peak_control, "ratio_QP")
mww_QP_table <- data.frame(
Statistic = QP_fire_stats$Statistic,
FIRE_QP = QP_fire_stats$Value,
CONTROL_QP = QP_ctrl_stats$Value
)
knitr::kable(
mww_QP_table,
digits = 4,
caption = "Mann–Whitney tests for Q_peak / P_peak: FIRE vs CONTROL"
)
Mann–Whitney tests for Q_peak / P_peak: FIRE vs
CONTROL
| Pre-fire median |
5.1589 |
260.0089 |
| Post-fire median |
12.5217 |
294.0789 |
| Median difference (post - pre) |
8.7056 |
36.1176 |
| 95% CI lower |
1.9247 |
-122.0362 |
| 95% CI upper |
25.0331 |
268.6274 |
| p-value |
0.0184 |
0.6171 |
## ============================================================
## 3) CHANGE-POINT TESTS (PETTITT) – FLASHINESS, Q/P, RR
## ============================================================
## 3a) Flashiness – FIRE site
pettitt_flash_stream <- run_pettitt(
wy_vec = flash_stream$wy,
value_vec = flash_stream$flash,
label = "Flashiness – FIRE site (stream)"
)
##
## Pettitt test for Flashiness – FIRE site (stream)
## Change-point WY: 2015
## K statistic : 26
## p-value : 0.3601846
kable(
data.frame(
Series = "Flashiness – FIRE",
Change_point_WY = pettitt_flash_stream$cp_wy,
p_value = pettitt_flash_stream$result$p.value
),
caption = "Pettitt Test – Flashiness (FIRE site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Flashiness (FIRE site)
|
Series
|
Change_point_WY
|
p_value
|
|
Flashiness – FIRE
|
2015
|
0.3601846
|
## 3b) Flashiness – CONTROL site
pettitt_flash_control <- run_pettitt(
wy_vec = flash_control$wy,
value_vec = flash_control$flash,
label = "Flashiness – CONTROL site"
)
##
## Pettitt test for Flashiness – CONTROL site
## Change-point WY: 2010 2011
## K statistic : 26
## p-value : 0.3601846
kable(
data.frame(
Series = "Flashiness – CONTROL",
Change_point_WY = pettitt_flash_control$cp_wy,
p_value = pettitt_flash_control$result$p.value
),
caption = "Pettitt Test – Flashiness (CONTROL site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Flashiness (CONTROL site)
|
Series
|
Change_point_WY
|
p_value
|
|
Flashiness – CONTROL
|
2010
|
0.3601846
|
|
Flashiness – CONTROL
|
2011
|
0.3601846
|
## 3c) Q_peak / P_peak – FIRE site
pettitt_ratio_stream <- run_pettitt(
wy_vec = peak_stream$wy,
value_vec = peak_stream$ratio_QP,
label = "Q_peak / P_peak – FIRE site (stream)"
)
##
## Pettitt test for Q_peak / P_peak – FIRE site (stream)
## Change-point WY: 2011
## K statistic : 34
## p-value : 0.1066308
kable(
data.frame(
Series = "Q_peak / P_peak – FIRE",
Change_point_WY = pettitt_ratio_stream$cp_wy,
p_value = pettitt_ratio_stream$result$p.value
),
caption = "Pettitt Test – Q/P Ratio (FIRE site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Q/P Ratio (FIRE site)
|
Series
|
Change_point_WY
|
p_value
|
|
Q_peak / P_peak – FIRE
|
2011
|
0.1066308
|
## 3d) Q_peak / P_peak – CONTROL site
pettitt_ratio_control <- run_pettitt(
wy_vec = peak_control$wy,
value_vec = peak_control$ratio_QP,
label = "Q_peak / P_peak – CONTROL site"
)
##
## Pettitt test for Q_peak / P_peak – CONTROL site
## Change-point WY: 2012
## K statistic : 16
## p-value : 1
kable(
data.frame(
Series = "Q_peak / P_peak – CONTROL",
Change_point_WY = pettitt_ratio_control$cp_wy,
p_value = pettitt_ratio_control$result$p.value
),
caption = "Pettitt Test – Q/P Ratio (CONTROL site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Q/P Ratio (CONTROL site)
|
Series
|
Change_point_WY
|
p_value
|
|
Q_peak / P_peak – CONTROL
|
2012
|
1
|
## 3e) Runoff Ratio – FIRE site
pettitt_RR_stream <- run_pettitt(
wy_vec = RR_stream$wy,
value_vec = RR_stream$RR,
label = "Runoff Ratio – FIRE site (stream)"
)
##
## Pettitt test for Runoff Ratio – FIRE site (stream)
## Change-point WY: 2012
## K statistic : 42
## p-value : 0.02281754
kable(
data.frame(
Series = "Runoff Ratio – FIRE",
Change_point_WY = pettitt_RR_stream$cp_wy,
p_value = pettitt_RR_stream$result$p.value
),
caption = "Pettitt Test – Runoff Ratio (FIRE site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Runoff Ratio (FIRE site)
|
Series
|
Change_point_WY
|
p_value
|
|
Runoff Ratio – FIRE
|
2012
|
0.0228175
|
## 3f) Runoff Ratio – CONTROL site
pettitt_RR_control <- run_pettitt(
wy_vec = RR_control$wy,
value_vec = RR_control$RR,
label = "Runoff Ratio – CONTROL site"
)
##
## Pettitt test for Runoff Ratio – CONTROL site
## Change-point WY: 2014
## K statistic : 22
## p-value : 0.5861141
kable(
data.frame(
Series = "Runoff Ratio – CONTROL",
Change_point_WY = pettitt_RR_control$cp_wy,
p_value = pettitt_RR_control$result$p.value
),
caption = "Pettitt Test – Runoff Ratio (CONTROL site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Runoff Ratio (CONTROL site)
|
Series
|
Change_point_WY
|
p_value
|
|
Runoff Ratio – CONTROL
|
2014
|
0.5861141
|
## 3g) Mean daily discharge – FIRE site (stream)
pettitt_meanQ_stream <- run_pettitt(
wy_vec = flash_stream$wy,
value_vec = flash_stream$meanQ,
label = "Mean Q – FIRE site (stream)"
)
##
## Pettitt test for Mean Q – FIRE site (stream)
## Change-point WY: 2012
## K statistic : 42
## p-value : 0.02281754
kable(
data.frame(
Series = "Mean Q – FIRE",
Change_point_WY = pettitt_meanQ_stream$cp_wy,
p_value = pettitt_meanQ_stream$result$p.value
),
caption = "Pettitt Test – Mean Daily Discharge (FIRE site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Mean Daily Discharge (FIRE site)
|
Series
|
Change_point_WY
|
p_value
|
|
Mean Q – FIRE
|
2012
|
0.0228175
|
## 3h) Mean daily discharge – CONTROL site
pettitt_meanQ_control <- run_pettitt(
wy_vec = flash_control$wy,
value_vec = flash_control$meanQ,
label = "Mean Q – CONTROL site"
)
##
## Pettitt test for Mean Q – CONTROL site
## Change-point WY: 2013 2014
## K statistic : 20
## p-value : 0.7252615
kable(
data.frame(
Series = "Mean Q – CONTROL",
Change_point_WY = pettitt_meanQ_control$cp_wy,
p_value = pettitt_meanQ_control$result$p.value
),
caption = "Pettitt Test – Mean Daily Discharge (CONTROL site)"
) |>
kable_styling(full_width = FALSE, bootstrap_options = "striped")
Pettitt Test – Mean Daily Discharge (CONTROL site)
|
Series
|
Change_point_WY
|
p_value
|
|
Mean Q – CONTROL
|
2013
|
0.7252615
|
|
Mean Q – CONTROL
|
2014
|
0.7252615
|