CCVA_Report.Rmd

---
# pulled some of this from Amber's climate report .Rmd for us to modify: https://github.com/nationalparkservice/CCRP-Exposure-Reports/blob/main/CCExposure.Rmd
params:
#  park_code: "BRCA" # eventually we will have a render script where we can map over park parameter names here
#  park_name: "Bryce Canyon National Park"
theme: united
header-includes:
   - \usepackage{booktabs}
output:
  word_document:
    reference_docx: template_CCRP.docx
    fig_width: 8
    fig_caption: true
    #toc: yes
    #toc_depth: 2
bibliography: references.bib
#csl: ecology.csl
editor_options: 
  markdown: 
    wrap: 79
  chunk_output_type: console
always_allow_html: true
---

```{r setup, include=FALSE}
# set chunk defaults
knitr::opts_chunk$set(echo = FALSE,
                      warning = FALSE,
                      error = FALSE,
                      message = FALSE,
                      fig.width = 7)
                      
                      
# read in libraries
source("setup.R")

# set parameters

# path to data folder (from project directory)
path <- "data/park/"

# park code for this report
park <- "BRCA"

# read in data
load(paste0(path, "/", park, "/", park, "_report_data.RData"))
## terra object issue when saving as .RData, so saved as separate .tif
wbm_30y <- terra::rast(paste0(path, "/", park, "/", park, "_wbm_30yr_annual_rasters.tif"))

# match Amber's 'Appendix_Script.R' to create env vars to use for report: https://github.com/nationalparkservice/CCRP-Exposure-Reports/blob/main/Appendix_Script.R
#source("Appendix-Script.R")
```

```{r nps_header, echo=FALSE, fig.width=7}
# CM: I for the life of me cannot get an image to fit to the width of the page and be centered
knitr::include_graphics("data/all/Report_Figures/NPS_Header.png")
```

# `r report_name`

::: {custom-style="Emphasis"}
NRSS Climate Change Response Program, NRSS Water Resources Division, Colorado
State University
:::

::: {custom-style="H1 box"}
Key Findings
:::

::: {custom-style="Text box"}
High level summary of key findings. Metrics such as projected change in runoff,
etc. -- need to write report first.
:::

**Recommended Citation**

```{r park_photo, echo=FALSE, out.width = 950}

knitr::include_graphics("data/all/Report_Figures/Park_Photos/BRCA.jpg")

```

\newpage

## Project Background

Rising temperatures, altered precipitation patterns, stronger storms, and other
climate changes are evident across America's national parks [@monahan2014climate]. 
These changes threaten the parks' natural, cultural, and built resources, which in turn can
impact opportunities to visit and recreate in these spaces. The National Park
Service (NPS) recommends that parks incorporate climate considerations into
planning for the future, emphasizing the need to develop forward-looking goals
that consider a range of divergent, plausible, and relevant future climatic
conditions.

This climate change vulnerability assessment (CCVA) focuses on how water
supplies at **`r park_name`** might be affected by projected shifts in climate
and hydrology. By selecting climate models that represent a range of possible
climate futures at `r park_name_short` from "Warm Wet" to "Hot Dry", the spread
of possible futures can be considered. Using the temperature, humidity, and
precipitation projections generated from these climate models, we then
calculate key hydrology-related variables that impact water supplies, including
snow, rain, soil moisture, evapotranspiration, and runoff. These projections
are intended to help parks make short- and long-term management decisions that
avoid surprises and costly mistakes.

## Regional Changes in Climate and Hydrology

The Intermountain West is already experiencing warming temperatures, shorter
winters with less snowfall, and increased drought. This is leading to more
extreme heat events, a higher risk of wildfires, and generally less water
availability. Climate-related changes in water may impact park water resources
in many ways, including:

-   **Reduced snowpack and earlier melt:** Rising temperatures will lead to a
    decrease in winter snowfall and an earlier spring melt. This will result in
    less water being stored in snowpack, which is a critical source for
    replenishing rivers and streams throughout the year.

-   **Increased droughts and water scarcity:** With less snowpack and altered
    precipitation patterns, droughts are expected to become more frequent and
    severe. This will lead to decreased streamflow, lower reservoir levels, and
    potential water shortages for agriculture, municipalities, and ecosystems.

-   **Shifts in Streamflow Timing:** Earlier snowmelt will cause peak
    streamflow to occur earlier in the spring, potentially leading to increased
    likelihood of flooding during that time. However, with less overall water
    stored as snow, streamflow will likely be lower during the summer months
    when water demands are typically highest.

-   **Changes in Precipitation Patterns:** While overall precipitation may not
    change dramatically, the pattern of precipitation is expected to shift.
    There might be an increase in extreme precipitation events (heavy rain and
    snowfall) interspersed with longer dry periods. This can lead to flash
    flooding and soil erosion during heavy rain events, while longer dry
    periods exacerbate drought conditions.

-   **Potential Impacts on Water Quality:** Lower streamflow can lead to
    increased water temperature and reduced dilution of pollutants. This can
    have negative consequences for aquatic ecosystems and water quality for
    human consumption.

-   **Feedback loops:** Complex feedbacks between climate and visitation,
    vegetation, wildfires, land management, and human water use, among others,
    have have potential to further impact water supplies in ways that are
    difficult or impossible to predict.

## `r park_name_short` Water Supply

```{r supply_text, echo = FALSE}

# create supply names text
a_sup <- active_supply_table$water_system_name2

if (length(a_sup) > 1) {
  a_sup <- paste(c(a_sup[1:length(a_sup)-1],paste0("and ",a_sup[length(a_sup)])), 
                      collapse = ", ", sep = "") 
}

num_2_word <- function(number) {
  # takes in an integer number between 1 and 20 and returns the word
  numbers_dict <- c(
    one = 1, two = 2, three = 3, four = 4, five = 5,
    six = 6, seven = 7, eight = 8, nine = 9, ten = 10,
    eleven = 11, twelve = 12, thirteen = 13, fourteen = 14,
    fifteen = 15, sixteen = 16, seventeen = 17, eighteen = 18,
    nineteen = 19, twenty = 20)
  return(names(numbers_dict[numbers_dict == number]))
}

if (nrow(active_supply_table) == 1) {
  sup_text <- paste0(a_sup," serves as the primary active water supply")
} else if (nrow(active_supply_table) > 1) {
  sup_text <- paste0(a_sup, " serve as the primary active water supplies")
}

sup_text2 <- rep(list(),nrow(active_supply_table))
for (i in 1:nrow(active_supply_table)) {
  sup_text2[i] <- paste0("**", a_sup[i] , ":** ", active_supply_table$description[i])
} 

source_types_at_park <- unique(active_supply_table$system_type_s)

gw_source <- "Groundwater is currently the primary source of water at the park. Groundwater is stored in underground aquifers that are recharged when surface water seeps into the ground.  Rain, streams, and other surface-water sources contribute to replenishing aquifers. To assess how climate change may impact this water supply, we must therefore consider both surface and groundwater systems."

sw_source <- "Surface water is currently the primary source of water at the park.  Surface water is sourced by rain and snow that travels at and beneath the surface towards drainages.  Many processes can intercept or diminish this water source, reducing the water that flows into streams.  To assess how climate change may impact this water supply, we must consider how the sources and interceptors of surface water might change."

sw_and_gw <- "Surface water and groundwater are currently the primary sources of water at the park.  This water is sourced by rain and snow that either flows towards streams or infiltrates into aquifers. Many processes can impact these flow paths, affecting the volume of water in each of these stores.  To assess how climate change may impact this water supply, we must consider how the many flow paths will change."

  if (grepl("Groundwater",paste(source_types_at_park, collapse = ""))) {
    source_text <- gw_source
  } else if (grepl("Surface Water",paste(source_types_at_park, collapse = ""))) {
    source_text <- sw_source
  } else if (grepl("Groundwater",paste(source_types_at_park, collapse = "")) &
        grepl("Surface Water",paste(source_types_at_park, collapse = ""))) {
    source_text <- sw_and_gw
  }
```

The National Park Service has `r num_2_word(nrow(POD_all))` points of diversion
(PODs) (**Figure 1**) and `r num_2_word(nrow(supply_table))` water supply
systems in or near `r park_name` (**Table 1**). `r sup_text` for
`r park_name_short`.

`r sup_text2`

`r source_text`

```{r park_map, fig.height=4}

tmap_mode(mode = c("plot"))
#tmap_mode(mode = c("view"))

map_area <- park_boundary %>%
  st_union(., nwis_stream) %>%
  st_union(., POD_all) %>%
  st_buffer(., 0.01) %>%
  st_bbox()

#tile_maps <- get_tiles(x = park_boundary, provider =  "Esri.WorldTerrain")
#tile_maps <- get_tiles(x = park_boundary, provider =  "Esri.NatGeoWorldMap")
#tile_maps <- get_tiles(x = park_boundary, provider =  "CartoDB.Positron")
#tile_maps <- get_tiles(x = park_boundary, provider =  "Esri.WorldTopoMap")
tile_maps <- get_tiles(x = map_area, provider =  "OpenStreetMap.HOT")
#tile_maps <- get_tiles(x = park_boundary, provider =  "OpenStreetMap")
#tile_maps <- get_tiles(x = park_boundary, provider =  "Esri.WorldImagery")
tile_maps2 <- get_tiles(x = map_area, provider =  "Esri.WorldShadedRelief")



## Create map
tm_shape(tile_maps, bbox = map_area) + 
  tm_rgb() +
tm_shape(tile_maps2, bbox = map_area) + 
  tm_rgb(alpha = 0.2) +
  tm_shape(park_boundary) +
  tm_polygons(col = "seagreen4", alpha = 0.25) +
  tm_shape(st_union(watersupply_watershed)) +
  tm_borders(col = "dodgerblue", lwd = 3) +
  tm_shape(POD_all) +
  tm_dots(col = "#BED558",
          alpha = 1,
          size = 0.4,
          shape = 21) +
  tm_shape(nwis_stream) +
  tm_symbols(col = "#8A5082",
          alpha = 1,
          size = 0.5,
          shape = 24) +
  tm_shape(nwis_groundwater) +
  tm_symbols(col = "#8A5082",
          alpha = 1,
          size = 0.5,
          shape = 22) +
  tm_shape(POD_supply, alpha = 0.25) +
  tm_dots(col = "dodgerblue", 
          border.col = "black", 
          border.lwd = 0.1,
          shape = 21,
          alpha = 1,
          size = 0.2) +
  tm_compass() +
  tm_scale_bar() +
  tm_add_legend("line", col = "dodgerblue", lwd = 3, labels = paste0(POD_supply$name[1]," Watershed")) +
  tm_add_legend("symbol", col = "dodgerblue",alpha = .9, labels = paste0(POD_supply$name[1]," Supply Point")) +
  tm_add_legend("symbol", col = "#BED558", alpha = 1, labels = "NPS-Owned PODs") +
  tm_add_legend("symbol", col = "#8A5082", alpha = 1, shape = 24, labels = "NWIS Stream Gauges") +
  tm_add_legend("symbol", col = "#8A5082", alpha = 1, shape = 21, labels = "NWIS Wells") +
  tm_layout(
      legend.text.size =.8,
      frame = FALSE,
      #legend.position = c(0.6, 0.20),
      legend.position = c("left","center"),
      #legend.text.size = .5,
      legend.outside = TRUE,
      #legend.outside.size = 0.5,
      outer.margins = 0,
      inner.margins = 0)



#tileURL  <- 'https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg&'

#aoaMap <- leaflet::leaflet() %>%
#  leaflet::addTiles(urlTemplate = tileURL) %>%
#  leaflet::addPolylines(data = park_boundary, label = "TEST", 
#                        color = "blue", weight = 3, opacity = 1)
#  htmltools::tagList(aoaMap %>%
#                     leaflet::addControl("TEST", position = "topright"))


```

***Figure 1**. Map water supply location, water supply watershed, and other
NPS-Owned points of diversion (PODs) at `r park_name`.\
*

***Table 1**. Summary table for listed water supply systems at
`r park_name_short`.*

```{r water_supply}

flextable(supply_table %>%
        dplyr::select(c("System Name" = water_system_name2, 
                        "System Type" = system_type_s,
                        "PWSS #" = pws_number,
                        "Classification" = classification),
                        "Active (y/n)" = active)) %>% 
bold(bold = TRUE, part = "header") %>% 
autofit() %>% 
fit_to_width(8)

```

\

## How We Assess Water Supplies: Water Budgets

Precipitation that falls to the ground as rain or snow can follow various
paths. Some infiltrates the soil, recharging groundwater and potentially
sustaining plant life. Some evaporates back into the atmosphere, especially
from surfaces like leaves and open water. A small portion might be stored as
surface water in lakes, ponds, or wetlands. The rest flows across the land
surface as runoff, eventually feeding streams, rivers, and ultimately the ocean
(**Figure 2**).

```{r wbm_schematic, echo=FALSE, fig.width = 4}
# CM: I for the life of me cannot get an image to fit to the width of the page and be centered
knitr::include_graphics("data/all/Report_Figures/NPS_WBM2.png")
```

***Figure 2**. Schematic illustration of the water budget.*

Each path represents a component of the "water budget". A water budget is like
a financial budget for water. It tracks all the water that enters (rain, snow),
exits (evapotranspiration, runoff), and gets stored (recharged groundwater,
soil moisture) in the system. The runoff, recharge, and storage components of
the water budget make up the vast majority park water supplies (i.e.,streams,
aquifers, and springs). By tracking these components, we can quantify changes
in water sources that feed park water supplies.

While analyzing historical data provides a valuable starting point for
estimating future water supplies, historical data has inherent limitations.
Climate change can introduce complex feedback loops that introduce new trends
and render historical trends unreliable. To address this challenge, Global
Climate Models (GCMs) are used along with water balance models to provide more
accurate projections for future climatic and hydrologic conditions. Because
many climate models exist, we use predetermined methods for selecting models
that are reliable and encompass the range for all model projections.

For `r park_name_short`, we select two climate models that represent "Warm Wet"
and "Hot Dry" futures that are plausible for the Intermountain Region.
Temperature, precipitation, and humidity from historic datasets and selected
models are then input to water balance models to generate past and future
estimates of park water budgets. Runoff is a key component of the water budget,
serving as a water source for stream discharge and groundwater recharge. Runoff
can be thought of as the rain and snow that isn't used by plants or stored in
the soil. **Figure 3** plots fractional components of the water budget for
historic and future scenarios at `r park_name_short`. Note how rain, snow,
evapotranspiration (AET) and runoff are predicted to change.

```{r WB_pie, echo = FALSE, fig.height = 2, fig.width = 5}
ggplot(annual_stats %>%
         mutate(runoff_in_mean = runoff_in_mean *365.25) %>%
         pivot_longer(cols = c(rain_in_mean, snow_in_mean, runoff_in_mean, 
                               aet_in_mean),#, soil_water_in_mean), 
                      values_to = "vals", names_to = "vars") %>%
         mutate(vars2= factor(vars, levels=c('rain_in_mean','snow_in_mean',
                                             'aet_in_mean','runoff_in_mean'))), 
       aes(x = "", y = vals, fill = vars2)) + 
  geom_bar(stat = "identity", width = 1, color = "white") + 
  coord_polar("y",start=0) +
  scale_fill_manual("",
              values = c("#0072B2","#56B4E9","#E69F00","#009E73","#CC79A7"),
                    labels = c("Rain", "Snow", "AET", "Runoff")) +
  facet_wrap(~cf, scales = "free_y") +
    theme_void() + 
  scale_y_discrete(expand = c(0, 0)) +
  theme(plot.margin = margin(0,.25,0,0,"cm"))
```

***Figure 3.** Water budget component fractions for the `r park_name_short`
water-supply watershed*

```{r wbm_text, echo = FALSE}

# Round magnitude for printing
magnitude <- magnitude %>%
  mutate(slope = round(slope,3),
         unit = case_when(grepl("_f",vars) ~ "degrees F",
                          grepl("_in",vars) ~ "inches",
                          grepl("days",vars) ~ "days",
                          grepl("TotalVisitors_all",vars) ~ "people",
                          grepl("S_D_ratio",vars) ~ ""))

mag_hist <- magnitude[magnitude$cf == "Historical",] %>%
  mutate(dir_text = 
           case_when(
             dir == "increases" ~ paste0("increased on average by ",round(slope,2)," ",unit," per year"),
             dir == "decreases" ~ paste0("decreased on average by ",round(slope,2)," ",unit," per year"),
             dir == "no significant changes" ~ "exhibited no significant changes"))

mag_hot_dry <- magnitude[magnitude$cf == "Hot Dry",] %>%
  mutate(dir_text = 
           case_when(
             dir == "increases" ~ paste0("increases on average by ",round(slope,2)," ",unit," per year."),
             dir == "decreases" ~ paste0("decreases on average by ",round(slope,2)," ",unit," per year."),
             dir == "no significant changes" ~ "exhibits no significant changes."))

mag_warm_wet <-magnitude[magnitude$cf == "Warm Wet",] %>%
  mutate(
    dir_text = 
      case_when(dir == "increases" ~ paste0("increases on average by ",round(slope,2)," ",unit," per year"),
                dir == "decreases" ~ paste0("decreases on average by ",round(slope,2)," ",unit," per year"),
                dir == "no significant changes" ~ "exhibits no significant changes"))

  

```

## Changes in climate and water budgets at `r park_name_short`

As global temperatures rise and precipitation patterns change, the delicate
balance of water within watersheds is being disrupted, leading to shifts in the
timing and magnitude of water budget components.

At `r park_name_short`, the historical period (1980-2010) recorded various
changes in the water budget. Temperature
`r mag_hist[mag_hist$vars == "tavg_f",]$dir_text`, rain
`r mag_hist[mag_hist$vars == "rain_in",]$dir_text`, and snow
`r mag_hist[mag_hist$vars == "snow_in",]$dir_text`. At the same time,
evapotranspiration `r mag_hist[mag_hist$vars == "aet_in" ,]$dir_text` and
runoff `r mag_hist[mag_hist$vars == "runoff_in" ,]$dir_text`. Additionally, the
peak spring runoff
`r spring_runoff_yearly_change[spring_runoff_yearly_change$cf == "Historical",]$dir`
by on average
`r spring_runoff_yearly_change[spring_runoff_yearly_change$cf == "Historical",]$slope_round`
days per year while peak snow cover has been
`r spring_swe_yearly_change[spring_swe_yearly_change$cf == "Historical",]$dir`
by on average
`r spring_swe_yearly_change[spring_swe_yearly_change$cf == "Historical",]$slope_round`
days per year.

Similarly, projected climate futures show additional changes to the water
budget. **"Warm Wet"** projections for the time period from 2025-2017 include
temperature `r mag_warm_wet[mag_warm_wet$vars == "tavg_f",]$dir_text`, rain
`r mag_warm_wet[mag_warm_wet$vars == "rain_in",]$dir_text`, snow
`r mag_warm_wet[mag_warm_wet$vars == "snow_in",]$dir_text`, evapotranspiration
`r mag_warm_wet[mag_warm_wet$vars == "aet_in",]$dir_text`, and runoff
`r mag_warm_wet[mag_warm_wet$vars == "runoff_in",]$dir_text`. Additionally,
peak spring runoff will be
`r spring_runoff_yearly_change[spring_runoff_yearly_change$cf == "Warm Wet",]$dir`
and peak snow cover will be
`r spring_swe_yearly_change[spring_swe_yearly_change$cf == "Warm Wet",]$dir`.

**"Hot Dry"** projections for the time period from 2025-2070 similarly include
temperature `r mag_hot_dry[mag_hot_dry$vars == "tavg_f",]$dir_text`, rain
`r mag_hot_dry[mag_hot_dry$vars == "rain_in",]$dir_text`, snow
`r mag_hot_dry[mag_hot_dry$vars == "snow_in",]$dir_text`, evapotranspiration
`r mag_hot_dry[mag_hot_dry$vars == "aet_in",]$dir_text`, and runoff
`r mag_hot_dry[mag_hot_dry$vars == "runoff_in",]$dir_text`. Additionally, peak
spring runoff will be
`r spring_runoff_yearly_change[spring_runoff_yearly_change$cf == "Hot Dry",]$dir`
and peak snow cover will be
`r spring_swe_yearly_change[spring_swe_yearly_change$cf == "Hot Dry",]$dir`.

```{r wbm_timeseries, fig.width=6, fig.height=3.1, echo = FALSE, message = FALSE, warning = FALSE}

ggplot(centroid_annual %>% ungroup() %>%
         dplyr::select(y, cf, tavg_f, snow_in, rain_in, runoff_in, aet_in) %>%
         pivot_longer(cols = -c(cf,y), values_to = "vals", names_to = "vars") %>%
         mutate(name = case_when(vars == "runoff_in" ~ "Runoff (in)",
                          vars == "aet_in" ~ "AET (in)",
                          vars == "rain_in" ~ "Rain (in)",
                          vars == "tavg_f" ~ "Temperature (F)",
                          vars == "snow_in" ~ "Snow (in)")), 
       aes(x = y, y = vals, color = cf)) + 
  geom_line(size = 0.5) + 
  stat_smooth (geom="line", alpha=0.5, size=.5, span=1) +
  #geom_smooth(span = .2, size = 1, se = FALSE, alpha = .1) +
  labs(x = "Year", y = "") +
  scale_color_manual("",values = c("black","red","dodgerblue3")) +
  theme_bw() +
  facet_wrap(~ factor(name, levels = c('Temperature (F)', 'Rain (in)', 
                                       'Snow (in)', 'AET (in)', 
                                       'Runoff (in)')), 
             scales = "free", ncol = 3) +  
   theme(legend.position=c(.8, .25),
         legend.title=element_blank(),
         #legend.direction="horizontal",
      strip.background=element_rect(colour="black",
                                    fill="ivory3"))
```

***Figure 4.** Changes to water budget components at `r park_name`.*

A healthy water budget relies on a predictable seasonal cycle. Changes in
precipitation or evaporation during specific times (e.g., earlier snowmelt,
increased summer evaporation) disrupt this cycle. This can lead to water
shortages, flooding risks, ecosystem impacts, and infrastructure challenges.
Even if the total annual water budget remains constant, changes during specific
seasons can have cascading effects on water availability, ecosystems, and
infrastructure. **Figure 5** shows changes between historic and projected water
budgets throughout the year, highlighting how different climate futures may
disrupt seasonal water cycles.

```{r wbm_doy, fig.width=6, fig.height=3.1, echo = FALSE, message = FALSE, warning = FALSE}
ggplot(doy %>%
         filter(vars %in% c("tavg_f", "aet_in", "snow_in", 
                            "rain_in", "runoff_in")), 
       aes(x = ym, y = vals, color = cf)) + 
  geom_line(size = 0.5) + 
  geom_smooth(span = .2, size = 1, se = FALSE) +
  scale_x_date(date_breaks = "2 month", date_labels =  "%b")  + 
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~ factor(name, levels = c('Temperature (F)', 'Rain (in)', 
                                      'Snow (in)', 'AET (in)', 'Runoff (in)')), 
             scales = "free", ncol = 3) +
  labs(x = "", y = "") +
  scale_color_manual("",values = c("black","red","dodgerblue3")) +
  theme(legend.position=c(.8, .25),
         legend.title=element_blank(),
         #legend.direction="horizontal",
      strip.background=element_rect(colour="black",
                                    fill="ivory3"))
```

***Figure 5.** Historical and projected shifts in water budget components by
day of year.*

Shifts in precipitation and snowmelt patterns combined with changes in
evapotranspiration can lead to more extreme events. For the water supply, we're
especially interested in how frequently extreme runoff events will occur, for
example, high flow events (floods) and low flow events (droughts). At
`r park_name_short`, historical low runoff events (\< 90th percentile) occurred
during
`r round(100*streaks_summary[streaks_summary$cf == "Historical",]$freq_low,1)`%
of years while high runoff events (\> 90th percentile) occurred during
`r 100*streaks_summary[streaks_summary$cf == "Historical",]$freq_high`% of
years. In the future, extreme runoff events are likely to increase (**Figure
6**). Projected runoff values for the "Warm Wet" scenario are predicted to
occur during
`r round(100*streaks_summary[streaks_summary$cf == "Warm Wet",]$freq_low,1)`%
of years and high runoff events during
`r round(100*streaks_summary[streaks_summary$cf == "Warm Wet",]$freq_high,1)`%
of years. For the "Hot Dry" scenario, low runoff events are predicted to occur
during
`r round(100*streaks_summary[streaks_summary$cf == "Hot Dry",]$freq_low,1)`% of
years and high runoff events during
`r round(100*streaks_summary[streaks_summary$cf == "Hot Dry",]$freq_high,1)`%
of years. Extreme events may also increase in frequency. For example,
historical high and low runoff events occurred at a maximum duration of
`r streaks_summary[streaks_summary$cf == "Historical",]$max_low` and
`r streaks_summary[streaks_summary$cf == "Historical",]$max_high` year(s)
respectively. In the "Warm Wet" future, they are projected to occur at a
maximum duration of
`r streaks_summary[streaks_summary$cf == "Warm Wet",]$max_low` and
`r streaks_summary[streaks_summary$cf == "Warm Wet",]$max_high` year(s) while
for the "Hot Dry" future, they are projected to occur at a maximum duration of
`r streaks_summary[streaks_summary$cf == "Hot Dry",]$max_low` and
`r streaks_summary[streaks_summary$cf == "Hot Dry",]$max_high` year(s).

```{r runoff_thresholds, fig.width=5, fig.height=3, echo = FALSE, message = FALSE, warning = FALSE}
a=ggplot() +
  geom_line(data = centroid_annual, aes(x = y, y = runoff_in, color = cf)) +
  geom_hline(yintercept = historic_probs$QLOW, linetype = "dashed") +
  geom_hline(yintercept = historic_probs$QHIGH, linetype = "dashed") +
  #geom_point(data = filter(streak_check, low_year == 1 | high_year == 1), 
  #           aes(x = y, y = runoff_in, color = cf)) +
  theme_minimal() +
  labs(x = "",y = "Runoff (in)",color = "") +
  scale_color_manual("",values = c("black","red","dodgerblue3"))

b = ggplot() +
  geom_point(data = streak_check, aes(x = y, y = low_year - 1, 
                                      color = cf, 
                                      fill = cf), 
             shape =22,cex = 2) +
  geom_point(data = streak_check, aes(x = y, y = high_year, 
                                      color = cf, fill = cf), shape = 22, cex = 2) +
  scale_y_continuous(breaks = c(0, 1),
    labels = c(expression(paste("< ", P[10], " Historic")), 
            expression(paste("> ", P[90], " Historic")))) +
  labs(x = "", y = "") +
  theme_minimal() +
  theme(legend.position = "none") + 
  scale_color_manual("",values = c("black","red","dodgerblue3")) + 
  scale_fill_manual("",values = c("black","red","dodgerblue3")) 

ggarrange(a, b, ncol =1, align = "v", heights = c(3,1))

```

***Figure 6.** Frequency of extreme runoff events at `r park_name`.*

# Water Demand

`r park_name_short`'s water demand encompasses all the water used within the
park, including water for drinking, sanitation, landscaping, ecosystems, and
even firefighting. While historic water demands are typically reported to state
agencies, future water demands must be estimated. This requires understanding
the drivers for park water use. Water use varies depending on park size, park
activities (e.g., construction projects), water system efficiency (e.g.,
leaks), and visitor use. Importantly, it can also vary due to climate change.
For example, rising temperatures may prompt each visitor to drink more water or
park visitation to require more irrigation. Climate change may also impact park
visitation. For example, parks with with milder shoulder seasons in the future
may see extended visitation periods, while those that experience more extreme
temperatures could face reduced visitor numbers.

Between the years 2000 and 2023, annual water use at `r park_name_short` varied
greatly from year to year
`r annual_use_summary[annual_use_summary$cf=="Historical",]$text`. During that
time, annual visitation
`r annual_visitor_summary[annual_visitor_summary$cf == "Historical", ]$text`
while per-person water use `r monthly_pp_demand$text` However, conservation
alone can never fully reduce water use to zero. We predict future water demands
using models based on climate and population changes, factors that are likely
to drive water use (**Figure 7**). In the "Warm Wet" scenario, projected
demands `r annual_use_summary[annual_use_summary$cf=="Warm Wet",]$text` while
in the "Hot Dry" scenario, they
`r annual_use_summary[annual_use_summary$cf=="Hot Dry",]$text`.

```{r demand, fig.width=5, fig.height=3, message = FALSE, warning = FALSE}

water_right_af_per_year <- POD_supply %>%
  distinct("WRNUM", .keep_all = TRUE) %>%
  pull(CFS) %>%
  sum(., na.rm = TRUE) * 86400 * 365 * 2.29569e-5
  
# Fix formatting

a <- ggplot(centroid_annual %>% drop_na(demand_pp_gal),aes(x = y, y = demand_pp_gal)) +
  stat_rollapplyr(align = "center", width = 10, FUN = mean, size = 1) +
  geom_point(alpha = 0.5) +
  labs(x = "", y = "Per-Person water use (af)") + 
  theme_bw()

# Fix formatting
b <- ggplot(centroid_annual %>% drop_na(TotalVisitors_all)) +
  geom_line(aes(x = y, y = TotalVisitors_all, color = cf), size = 1) +
  geom_ribbon(aes(x = y, ymin =  TotalVisitors_mlr_lwr, ymax =  TotalVisitors_mlr_upr, 
              color = cf, fill = cf), alpha = 0.2, linewidth = 0) + 
  #geom_line(aes(x = y, 
   #             y = TotalVisitors),
  #              color = "black") + 
  scale_color_manual("", values = c("black","red","dodgerblue3")) +
  scale_fill_manual("", values = c("white","red","dodgerblue3")) +
  labs(x = "", y = "Annual Park Visitors") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6)) +
  xlim(1980,2070)  + theme_bw()

c <- ggplot(centroid_annual %>% drop_na(demand_all)) +
  geom_line(aes(x = y, y = demand_all, color = cf), size = 1) +

geom_ribbon(aes(x = y, ymin =  MLR_demandA_lwr, ymax =  MLR_demandA_upr, 
              color = cf, fill = cf), alpha = 0.2, linewidth = 0) + 
  geom_line(aes(x = y, 
                y = use_acre_feet),
                color = "black") + 
  geom_hline(yintercept=water_right_af_per_year, linetype = "dashed", color = "gray40") +
  geom_text(aes(2000,
                water_right_af_per_year,
                label = "Annual Water Right", 
                vjust = -.5),check_overlap = TRUE,
            size = 3, color = "gray40") +
  scale_color_manual("", values = c("black","red","dodgerblue3")) +
  scale_fill_manual("", values = c("black","red","dodgerblue3")) +
  labs(x = "", y = "Annual Water Demand (af)") + 
  xlim(1980,2070) + theme_bw()#+ ylim(-100,350) + theme_bw()

ggarrange(a,b,c, ncol = 1, scales="free", align = "hv", common.legend = TRUE)

  
ggplot(centroid_annual) +
  geom_line(aes(x = y, y = demand_all, color = cf)) +
```

***Figure 7**. Historical and projected annual water demand for `r park_name`*

# Vulnerability

Parks face water supply vulnerability when the balance between available water
(supply) and water use (demand) is disrupted. A resilient system maintains a
steady supply that can meet current and future water needs. Challenges can
arise when either water supply diminishes or water demand rises. Even minor
changes in water supply can be overwhelmed by a sharp rise in water use,
leading to a vulnerable system.

To assess water-supply vulnerability at `r park_name_short`, we calculate the
supply-to-demand ratio (S/D) for historical and projected supply and demand
at the park. The S/D ratio provides a snapshot of potential stress on the
system, and is a direct indicator for how system vulnerability is trending
through time. A shrinking ratio, where demand approaches or surpasses supply,
suggests growing vulnerability that can signal potential water shortages,
increased competition for resources, or an urgent need to implement water
management strategies that ensure sustainable water use. Monitoring this ratio
over time can therefore help us anticipate and address potential water crises
before they occur.

At `r park_name_short`, the average S/D ratio for recent years (2000-2023) has
been `r round(annual_SD_summary[annual_SD_summary$cf == "Historical",]$mean,2)`
with a Sen's slope indicating that value is
`r annual_SD_summary[annual_SD_summary$cf == "Historical",]$text`. The "Hot
Dry" projections suggest a future S/D ratio that is
`r annual_SD_summary[annual_SD_summary$cf == "Hot Dry",]$text` (approximately
`r round(annual_SD_summary[annual_SD_summary$cf == "Hot Dry",]$percent_chg,1)`
% per year).\
The "Warm Wet" projections suggest an S/D ratio that is
`r annual_SD_summary[annual_SD_summary$cf == "Warm Wet",]$text` (approximately
`r round(annual_SD_summary[annual_SD_summary$cf == "Warm Wet",]$percent_chg,1)`
% per year.

```{r s_d_ratio, fig.width=5, fig.height=3, message = FALSE, warning = FALSE}


ggplot(centroid_annual, aes(x = y, y = S_D_ratio, color = cf)) +
  geom_line() +
  geom_smooth(method = "lm",aes(fill= cf)) +
  #geom_ribbon(aes(x = y, ymin =  S_D_ratio_lwr, ymax =  S_D_ratio_upr, 
  #            color = cf, fill = cf), alpha = 0.2, linewidth = 0) + 
  scale_color_manual("", values = c("black","red","dodgerblue3")) +
  scale_fill_manual("", values = c("black","red","dodgerblue3")) +
  #scale_fill_manual("", values = c("black","darkred","dodgerblue4")) +
  labs(x = "", y = "S/D Ratio") + 
  theme_bw()

```

***Figure 8**. Historical and projected supply-to-demand ratio (S/D ratio) for
`r park_name`.*

# Uncertainty

### Water rights limitations

In the western US, water rights are a looming challenge for future water
stress. Many water rights were historically allocated during times of greater
abundance and may not be flexible enough to adapt to a future marked by reduced
precipitation and a growing population. This inflexibility could exacerbate
water stress, creating conflicts between established rights.

### Spatial variability

Climate-related changes in water budgets will not be evenly distributed.
Understanding spatial variations is crucial for developing targeted water
management strategies and adapting to the challenges posed by climate change in
different regions. Things that can affect water budget response to climate
change include latitude, global wind and ocean patterns, local topography, and
landcover changes. In general, we can expect increased variability, greater
impacts in vulnerable regions, and an uneven distribution of impacts. **Figure
9** shows how runoff is projected to change throughout the area surrounding
`r park_name_short`.

```{r runoff_change, fig.width=7, fig.height = 3, echo = FALSE, message = FALSE, warning = FALSE}
# create change maps, CM single example: 2055  hot dry change for runoff

# calculate change rasters (%) for runoff
runoff_change_rasters <- map2(rep(scen, 2), c(rep("2040_2069", 2), rep("2070_2099", 2)), function(x, y) {
  (wbm_30y[[grepl(pattern = paste(select_cfs[select_cfs$CF == x, GCM], "runoff", y, sep = "_"),
                  x = names(wbm_30y))]] - wbm_30y[[grepl(
                    pattern = paste("gridmet_historical", "runoff", "1981_2010", sep = "_"),
                    x = names(wbm_30y)
                  )]]) / wbm_30y[[grepl(
                    pattern = paste("gridmet_historical", "runoff", "1981_2010", sep = "_"),
                    x = names(wbm_30y)
                  )]] * 100
}) %>% terra::rast() %>% 
  crop(st_buffer(st_transform(park_boundary, crs(.)), 10000))

# create all figures
runoff_change_maps <- map(scen, function(x) {
  # map2(rep(scen, 2), c(rep("2040_2069", 2), 
  #                                          rep("2070_2099", 2)), function(x, y) {
  tm_shape(runoff_change_rasters[[grepl(
    pattern = paste(select_cfs[select_cfs$CF == x, GCM], "runoff", "2040_2069", sep = "_"),
    x = names(runoff_change_rasters)
  )]]) +
    tm_raster(
      style = "cont",
      # breaks = c(-1000, -100, 0, 100, 1000),
      # palette = c("red","red", "white", "dodgerblue","dodgerblue"),
      breaks = c(-100, -50, 0, 50, 100),
      palette = c("red", "white", "dodgerblue"),
      midpoint = 0,
      legend.reverse = TRUE,
      n = 6,
      legend.hist = TRUE,
      title = "",#paste0("Predicted Change in\nRunoff: ", y),
      legend.format = list(
        fun = function(z) {
          ifelse(z < -100, "", 
                 ifelse(z > 100, "", 
                        ifelse(z == -100, "≤ -100 %", 
                               ifelse(z == 100, "≥ 100 %",
                                      paste0(formatC(
                                        z, digits = 0, format = "f"
                                      ), " %")))))
          
        }
        
      )) +
    tm_shape(park_boundary) +
    tm_borders(col = "black", lwd = 2) +
    tm_shape(st_union(watersupply_watershed)) +
    tm_borders(col = "gold1", lwd = 3) +
    tm_scale_bar(position = c("center", "bottom"), lwd = 2, text.size = 1) +
    tm_add_legend("fill",
                  col = "gold1",
                  labels = "Water Supply Watershed")+
    tm_layout(
      title = paste0("Δ Runoff, ", x, "\n2040-2070"),
      title.size = 0.9,
      frame = FALSE,
      legend.outside = TRUE,
      legend.outside.size = 0.35,
      outer.margins = c(0.001, 0.001, 0.001, 0.001),  # Reduce outer margins
      inner.margins = c(0.001, 0.001, 0.001, 0.001),
      legend.position = c("left", "center"),  # Position legend to avoid title overlap
      #legend.outside.position = "bottom" # Reduce inner margins
      
    )
})

tmap_arrange(runoff_change_maps, ncol = 2, outer.margins = -0.003)
```

***Figure 9**. Map of projected change in runoff for both CFs and both future
time periods.*

# Recommendations

**Incorporate climate change projections:** At `r park_name_short`, the current
water supply watershed is expected Integrate climate change projections into
water resource planning to prepare for future scenarios of reduced
precipitation or increased evaporation.

**Implement water conservation measures:** Water conservation measures can be
implemented by park visitors and staff. Encourage everyone in the park to adopt
water-saving practices (low-flow fixtures, shorter showers, etc.) through
education and outreach programs. Reduce irrigated areas by planting native
vegetation. Assess the feasibility of rainwater harvesting and gray water reuse
for lawn and native plant irrigation.

**Upgrade infrastructure:** Many parks have aging water supply infrastructure
that can lead to wasted water and inefficient water use. Investing in leak
detection and repair programs can help minimize water losses and improve the
efficiency of delivery systems.

**Improve data collection:** Projections are only as good as the data that goes
into developing them. Improve metering and reporting of water use in the park.
Enhance monitoring of water sources (snowpack, streamflow, groundwater levels)
to gain a better understanding of water availability and trends.

\newpage

# Appendix

Sen's slope is a non-parametric statistic used to assess trends in time series
data. It represents the median change between pairs of data points over time. A
positive Sen's slope indicates an upward trend, and a negative Sen's slope
indicates a downward trend. The value of the slope signifies the magnitude of
the change. The following table of Sen's slopes provides an indication of trends
for various climate and water balance variables for the historical period of 
record and projected scenarios at `r park_name`.

***Table S1.** Sen's slope estimator for climate and hydrology variables.*

```{r Stats_Table}

flextable(magnitude %>%
     mutate(vars = gsub("_in$", "",vars) %>%
            sapply(., function(x) {
              paste0(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)), sep = " ")
              })) %>%
       mutate(vars = gsub("_", " ", vars)) %>%
       mutate(vars = gsub("gt",">",vars)) %>%
       mutate(vars = gsub("lt","<",vars)) %>%
       mutate(vars = gsub("swe","SWE",vars)) %>%
       mutate(vars = gsub("pcp", "precip",vars)) %>%
       mutate(vars = gsub("roff"," runoff", vars)) %>%
       mutate(vars = gsub("precip", " precip", vars)) %>%
       mutate(vars = gsub("95","P_95", vars)) %>%
       mutate(vars = gsub("05","P_05", vars)) %>%
       mutate(vars = gsub(" w "," with ",vars)) %>%
        group_by(vars) %>%
        dplyr::select(cf,"Variable" = vars,slope) %>%
        pivot_wider(names_from = cf, values_from = c(slope)))%>% 
colformat_num(j = c("Historical", "Hot Dry",
                             "Warm Wet"), 
i = 2, digits = 0) %>% 
  bold(bold = TRUE, part = "header") %>% 
  autofit() %>% 
  fit_to_width(8)
```

### Citations

1.  <https://www.nps.gov/subjects/climatechange/planning.htm>

2.  Lawrence, D. J., Runyon, A. N., Gross, J. E., Schuurman, G. W., &
    Miller, B. W. (2021). Divergent, plausible, and relevant climate futures
    for near-and long-term resource planning. Climatic Change, 167(3), 38.

3.  Tercek, M. T., Thoma, D., Gross, J. E., Sherrill, K., Kagone, S., &
    Senay, G. (2021). Historical changes in plant water use and need in the
    continental United States. Plos one, 16(9), e0256586.