Climatic Analysis of Telescope Peak Bristlecone Pine Mortality
Table of Contents
- 1. Introduction
- 2. Data
- 3. Historical Water Balance Averages
- 4. Simple Water Balance Model (Excel) (Penman)
- 5. Lutz et al. 2010 Water Balance Method
- 6. Hypotheses
- 7. Questions
- 8. Next Steps
- 9. Discussion
- 10. Bibliography
1. Introduction
This project extends the climate analysis conducted on the PILO mortality cite documented by (Bentz et al. 2022). We attempt to recreate their analysis here, as well as run an additional analysis following David Thoma's simple water balance model which uses daily Daymet timeseries and incorporates additional physical site parameters and therefore may provided a more advanced, and hopefully accurate, analysis of the site.
Our aim is to examine the mortality event occuring at Panamint Range (Telescope Peak), which was attributed to Mountain Pine Beetle (Bentz et al. 2022). Mortality events also occured at Wah Wah Mountains (Pinyon Ips), Silver Peak (unattributed) , and Hot Creek Range/Rawhide (unattributed). While the mortality event was not recorded until 2022, it is believed that it occured around 2017. Points from Snake Range and Mount Moriah are included as a location where PIFL mortality was noted without associated PILO mortality.
Four analyses were performed:
- 3 Use historical water balance summary layers from Mike Tercek
- Simple Water Balance Model (Excel) (Penman)
shareable_Excel_xlsm_water_balance_model_v3.xlsmSimple Water Balance model implemented in Excel provided by David Toma, using Daymet daily time series - Lutz et al. 2010 Water Balance Method Lutz 2010 Water Balance Model implemented in R following methods in (Bentz et al. 2022), which followed (Lutz, van Wagtendonk, and Franklin 2010), using PRISM 30yr normals
1.1. Map
2. Data
2.1. DEM
SRTM 30m DEM used to generate slope, aspect, elevation in QGIS 3.10. Some elevations were provided in the Bentz data set, existing values were used instead of generated ones when available.
2.2. Soil
SSURGO database
2.3. t50
merged_jennings.tif
3. Historical Water Balance Averages
Using summary layers generated by Mike Tercek http://screenedcleanedsummaries.s3-website-us-west-2.amazonaws.co/.
3.1. AET x CWD
library(tidyverse) library(here) library(gghighlight) mt_list <- c("PANAMINT RANGE", "WAH WAH MOUNTAINS", "SILVER PEAK RANGE", "HOT CREEK RANGE", "SNAKE RANGE", "SNAKE RANGE(MORIAH)") "HOT CREEK RANGE") # sampled in QGIS tbl <- read_csv("./sites_aet_cwd_2000_2019_historical_avg.csv") ggplot(data = tbl) + geom_point(mapping = aes(x = CWD1, y = AET1, color = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + ggtitle("AET x CWD Annual gridMET Historical 2000-2019")
4. Simple Water Balance Model (Excel) (Penman)
Data from David Thoma's Simple Water Balance model xlsm, processed using provided batch script. Daymet data obtained using extractdaymetbatch.py script. PET Values calculated with Penman method
4.1. Libraries
library(tidyverse) library(ggplot2) library(readxl) library(here) library(gghighlight) library(ggpubr)
4.2. Data Import and Cleanup
setwd("penman/annual") tbl <- list.files(pattern = "*.csv") %>% map_df(~read_csv(., n_max = 41, col_types = cols())) %>% # only read 41 rows (For annual reports) because csvs contain nonsynctactic row of averages at bottom, which creates NAs in df mutate(year = strtoi(`Row Labels`)) %>% group_by(site) setwd(here()) sites_tbl <- read_excel("sites.xlsx") mt_list <- c("PANAMINT RANGE", "WAH WAH MOUNTAINS", "SILVER PEAK RANGE", "HOT CREEK RANGE", "SNAKE RANGE", "SNAKE RANGE(MORIAH)")
4.3. Summary Table
summary_tbl <- tbl %>% summarise( mean_t = mean(`Average of T`), mean_P = mean(`Sum of P`), mean_AET= mean(`Sum of AET`), mean_D= mean(`Sum of D`), mean_GDD= mean(`Sum of GDD`), mean_SOIL= mean(`Average of SOIL`), ### mean_t_pre_2010 = mean(`Average of T`[year <= 2010]), mean_P_pre_2010 = mean(`Sum of P`[year <= 2010]), mean_AET_pre_2010 = mean(`Sum of AET`[year <= 2010]), mean_D_pre_2010 = mean(`Sum of D`[year <= 2010]), mean_GDD_pre_2010 = mean(`Sum of GDD`[year <= 2010]), mean_SOIL_pre_2010 = mean(`Average of SOIL`[year <= 2010]), ### mean_t_post_2010 = mean(`Average of T`[year > 2010]), mean_P_post_2010 = mean(`Sum of P`[year > 2010]), mean_AET_post_2010 = mean(`Sum of AET`[year > 2010]), mean_D_post_2010 = mean(`Sum of D`[year > 2010]), mean_GDD_post_2010 = mean(`Sum of GDD`[year > 2010]), mean_SOIL_post_2010 = mean(`Average of SOIL`[year > 2010]), ### Delta values pre and post 2010 dt = mean_t_post_2010 - mean_t_pre_2010, dP = mean_P_post_2010 - mean_P_pre_2010, dAET = mean_AET_post_2010 - mean_AET_pre_2010, dD = mean_D_post_2010 - mean_D_pre_2010, dGDD = mean_GDD_post_2010 - mean_GDD_pre_2010, dSOIL = mean_SOIL_post_2010 - mean_SOIL_pre_2010, ) %>% left_join(sites_tbl, by = c("site" = "Site")) %>% group_by(MT_range)
4.4. Figures
4.4.1. AET x CWD
summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_D, y = mean_AET, col = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + ggtitle("AET x CWD Penman Daymet Excel")
4.4.2. Pre 2010 Plot
cwd_plt_pre_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_pre_2010, y = mean_D_pre_2010, color = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + labs(x = "Mean Annual T (1980-2010)", y = "CWD") p_plt_pre_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_pre_2010, y = mean_P_pre_2010, color = MT_range)) + labs(x = "Mean Annual T (1980-2010)", y = "PPT") + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) aet_plt_pre_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_pre_2010, y = mean_AET_pre_2010, color = MT_range)) + labs(x = "Mean Annual T", y = "AET") + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) ggarrange(cwd_plt_pre_2010, p_plt_pre_2010, aet_plt_pre_2010, common.legend = TRUE) %>% annotate_figure(top = text_grob("Climate Space 1980-2010 Penman Daymet Excel"))
4.4.3. Post 2010 Plot
cwd_plt_post_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_post_2010, y = mean_D_post_2010, color = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + labs(x = "Mean Annual T (2011-2020)", y = "CWD") p_plt_post_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_post_2010, y = mean_P_post_2010, color = MT_range)) + labs(x = "Mean Annual T (2011-2020)", y = "PPT") + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) aet_plt_post_2010 <- summary_tbl %>% ggplot() + geom_point(mapping = aes(x = mean_t_post_2010, y = mean_AET_post_2010, color = MT_range)) + labs(x = "Mean Annual T (2011-2020)", y = "AET") + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) ggarrange(cwd_plt_post_2010, p_plt_post_2010, aet_plt_post_2010, common.legend = TRUE) %>% annotate_figure(top = text_grob("Climate Space 2011-2020 Penman Daymet Excel"))
4.4.4. Change Plot
summary_tbl %>% ggplot(aes(x = mean_D_pre_2010, y = mean_AET_pre_2010, color = MT_range)) + geom_arrow(aes(dx = dD, dy = dAET)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + labs(title = "Change in mean AET and CWD over periods from 2010-2021 and 1980-2010 Penman Daymet Excel", x = "CWD (mm)", y = "AET (mm)") + scale_mag()
4.5. Timeseries
4.5.1. Load Monthly Data
setwd("penman/monthly") monthly_tbl <- list.files(pattern = "*.csv") %>% map_df(~read_csv(., col_types = cols())) %>% separate(`Row Labels`, c("year", "month")) %>% mutate(date = as.Date(paste(year, month, "01", sep = "/"), formate="%Y/%m/%d")) %>% left_join(sites_tbl, by = c("site" = "Site")) %>% arrange(MT_range, date) %>% group_by(MT_range) setwd(here())
4.5.2. Summary Table
monthly_summary_tbl <- monthly_tbl %>% group_by(MT_range, year, month) %>% summarize(AET = mean(`Sum of AET`), PET = mean(`Sum of PET`), D = mean(`Sum of D`), SNOW = mean(`Sum of SNOW`), PACK = mean(`Max of PACK`), SOIL = mean(`Average of SOIL`), GDD = mean(`Sum of GDD`), T = mean(`Average of T`), P = mean(`Sum of P`), RAIN = mean(`Sum of RAIN`), date = date) %>% distinct()
4.5.3. AET
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = AET, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("AET Monthly Time Series Penman Daymet Excel")
4.5.4. CWD
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = D, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("CWD Monthly Time Series Penman Daymet Excel")
4.5.5. T
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = T, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("Temperature Monthly Time Series Penman Daymet Excel")
4.5.6. PPT
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = P, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("PPT Monthly Time Series Penman Daymet Excel")
4.5.7. RAIN
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = RAIN, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("RAIN Monthly Time Series Penman Daymet Excel")
4.5.8. SNOW
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = SNOW, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("SNOW Monthly Time Series Penman Daymet Excel")
4.5.9. PACK
monthly_summary_tbl %>% filter(year >= 2008, year <= 2018) %>% ggplot(mapping = aes(x = date, y = PACK, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("PACK Monthly Time Series Penman Daymet Excel")
4.5.10. SOIL
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = SOIL, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("SOIL Monthly Time Series Penman Daymet Excel")
4.5.11. GDD
monthly_summary_tbl %>% filter(year >= 2012) %>% ggplot(mapping = aes(x = date, y = GDD, color = MT_range)) + geom_line() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + theme(legend.position = "bottom") + scale_x_date(date_breaks = "1 year", date_labels = "%Y", date_minor_breaks = "1 month") + ggtitle("GDD Monthly Time Series Penman Daymet Excel")
4.6. Animations
library(gganimate) mortality_tbl <- tbl %>% left_join(sites_tbl, by = c("site" = "Site")) %>% ungroup() %>% mutate(mortality = if_else(MT_range %in% c("PANAMINT RANGE", "WAH WAH MOUNTAINS", "HOT CREEK RANGE", "SILVER PEAK RANGE"), TRUE, FALSE))
4.6.1. AET x CWD over time
p <- tbl %>% left_join(sites_tbl, by = c("site" = "Site")) %>% ggplot() + geom_point(mapping = aes(x = `Sum of D`, y = `Sum of AET`, col = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE, calculate_per_facet = TRUE) + labs(title = 'Year: {frame_time}') + transition_time(year) + ease_aes('linear') animate(p, renderer = ffmpeg_renderer()) anim_save("aet_x_cwd.mp4", animation = last_animation(), path = "img")
4.6.2. All Climate variables over time
p <- mortality_tbl %>% group_by(mortality) %>% arrange(year) %>% pivot_longer(c("ASPECT_QGIS", "SLOPE_QGIS", "Average of SOIL", "Average of T", "Elev_m", "Max of PACK", "Sum of AET", "Sum of D", "Sum of GDD", "Sum of MELT", "Sum of P", "Sum of PET", "Sum of RAIN", "Sum of SNOW", "Sum of W", "Sum of W - PET", "Sum of W-ET-/\\SOIL")) %>% transform(name = factor(name, levels = c("ASPECT_QGIS", "SLOPE_QGIS", "Average of SOIL", "Average of T", "Elev_m", "Max of PACK", "Sum of AET", "Sum of D", "Sum of GDD", "Sum of MELT", "Sum of P", "Sum of PET", "Sum of RAIN", "Sum of SNOW", "Sum of W", "Sum of W - PET", "Sum of W-ET-/\\SOIL"))) %>% ggplot() + geom_boxplot(mapping = aes(x = mortality, y = value)) + facet_wrap(~ name, scales = "free") + ## transition_states( ## year, ## transition_length = 2, ## state_length = 1 ## ) + labs(title = 'Year: {frame_time}') + transition_time(year) + ease_aes('linear') animate(p, renderer = ffmpeg_renderer()) anim_save("boxplot_time.mp4", animation = last_animation(), path = "img")
4.7. Modelling
## log regression attempt ## for (y in c(2000:2017)) { ## mortality_tbl <- tbl %>% ## ungroup() %>% ## left_join(sites_tbl, by = c("site" = "Site")) %>% ## filter(year == y) %>% ## mutate(mortality = if_else(MT_range == "PANAMINT RANGE", ## TRUE, ## FALSE)) %>% ## ungroup() %>% ## select(ASPECT_QGIS, ## SLOPE_QGIS, ## `Average of SOIL`, ## `Average of T`, ## Elev_m, ## `Max of PACK`, ## `Sum of AET`, ## `Sum of D`, ## `Sum of GDD`, ## `Sum of MELT`, ## `Sum of P`, ## `Sum of PET`, ## `Sum of RAIN`, ## `Sum of SNOW`, ## `Sum of W`, ## `Sum of W - PET`, ## `Sum of W-ET-/\\SOIL`, ## mortality) ## mortality_tbl$mortality <- factor(mortality_tbl$mortality) ## mylogit <- glm(mortality ~ ., data = mortality_tbl, family = "binomial"(link='logit')) ## print(summary(mylogit)) ## } ## mortality_tbl <- tbl %>% ## ungroup() %>% ## left_join(sites_tbl, by = c("site" = "Site")) %>% ## filter(year == ) %>% ## mutate(mortality = if_else(MT_range == "PANAMINT RANGE", ## TRUE, ## FALSE)) %>% ## ungroup() %>% ## select(ASPECT_QGIS, ## SLOPE_QGIS, ## `Average of SOIL`, ## `Average of T`, ## Elev_m, ## `Max of PACK`, ## `Sum of AET`, ## `Sum of D`, ## `Sum of GDD`, ## `Sum of MELT`, ## `Sum of P`, ## `Sum of PET`, ## `Sum of RAIN`, ## `Sum of SNOW`, ## `Sum of W`, ## `Sum of W - PET`, ## `Sum of W-ET-/\\SOIL`, ## mortality) ## mortality_tbl$mortality <- factor(mortality_tbl$mortality) ## mylogit <- glm(mortality ~ ., data = mortality_tbl, family = "binomial"(link='logit')) ## print(summary(mylogit))
5. Lutz et al. 2010 Water Balance Method
Following (Bentz et al. 2022), Implement water balance model as in (Lutz, van Wagtendonk, and Franklin 2010) and supplemental. Use 800m PRISM 30 yr normals.
5.1. Libraries
library(raster) library(ggplot2) library(tidyverse) library(lubridate) library(gghighlight) library(prism)
Be sure to set the download folder using prism_set_dl_dir().
prism_set_dl_dir("~/prismtmp")
5.2. Water Balance Functions
Implement water balance functions from (Lutz, van Wagtendonk, and Franklin 2010) supplemental
get_f <- function (tmean) { f <- case_when( tmean <= 0 ~ 0, tmean > 0 & tmean < 6 ~ 0.167 * tmean, tmean >= 6 ~ 1) return(f) } get_rain <- function (ppt, F) { return(F * ppt) } get_snow <- function (ppt, F) { return( (1 - F) * ppt ) } get_pack <- function (ppt, F, sp.0=NULL) { snowpack <- vector() sp.0 <- ifelse(!is.null(sp.0), sp.0, 0) for (i in 1:length(ppt)) { if (i == 1) { snowpack[i] = (1 - F[i])**2 * ppt[i] + (1 - F[i]) * sp.0 } else { snowpack[i] = (1 - F[i])**2 * ppt[i] + (1 - F[i]) * snowpack[i - 1] } } return(snowpack) } get_melt <- function (snow, pack, F, sp.0=NULL) { sp.0 <- ifelse(!is.null(sp.0), sp.0, 0) melt <- vector() for (i in 1:length(snow)) { if ( i == 1 ) { melt[i] = F[i] * (snow[i] + sp.0) } else { melt[i] = F[i] * (snow[i] + pack[i-1]) } } return(melt) } get_dl <- function (mon, days, Lat) { ## Get Daylength for all days in vector of months date <- paste("1980-", mon, "-", days, sep = "") yd <- yday(date) theta <- 0.2163108+2*atan(0.9671396*tan(0.00860*(yd-186))) P <- asin(0.39795 * cos(theta)) dl <- 24 - (24/pi) * acos((sin((0.8333 * pi)/180) + sin((Lat * pi) / 180) * sin(P))/(cos((Lat*pi)/180)*cos(P))) return(dl) } get_hl <- function (Lat, slope, aspect_f) { ## calculate heat load index multiplier Lat.rad <- (pi/180) * Lat slope.rad <- (pi/180) * slope HL <- 0.339 + 0.808 * (cos(Lat.rad) * cos(slope.rad)) - 0.196 * (sin(Lat.rad) * sin(slope.rad)) - 0.482 * (cos(aspect_f) * sin(slope.rad)) return(HL) } get_soil <- function (soil_max, w, pet, s.0=NULL) { s.0 = ifelse(!is.null(s.0), s.0, 0) soil <- vector() for (i in 1:length(pet)) { if ( i == 1 ) { soil[i] = pmin(soil_max[i], if (w[i] > pet[i]) { (w[i] - pet[i]) + s.0 } else { s.0 * (1 - exp(-(pet[i]-w[i])/soil_max[i])) }) } else { soil[i] = pmin(soil_max[i], if (w[i] > pet[i]) { (w[i] - pet[i]) + soil[i-1] } else { soil[i-1] * (1 - exp(-(pet[i]-w[i])/soil_max[i])) }) } } return(soil) } get_d_soil <- function (soil, s.0=NULL) { s.0 = ifelse(!is.null(s.0), s.0, 0) d_soil = soil - lag(soil, default = s.0) return(d_soil) } get_aet <- function (pet, d_soil, w) { aet <- vector() for (i in 1:length(pet)) { a <- min(pet[i], d_soil[i] + w[i]) aet[i] = if_else(a > 0, a, 0) } return(aet) }
5.3. Download PRISM Normals
Only needs to be run once after which data is saved to prism dl dir
get_prism_normals("ppt", "800m", annual = TRUE, keepZip = FALSE) get_prism_normals("ppt", "800m", mon = 1:12, keepZip = FALSE) get_prism_normals("tmean", "800m", annual = TRUE, keepZip = FALSE) get_prism_normals("tmean", "800m", mon = 1:12, keepZip = FALSE)
5.4. Data Import and Cleanup
points <- read_csv("./sites.csv") points.spdf <- SpatialPointsDataFrame(coords = points[,c('Lon', 'Lat')], data = points, proj4string = CRS("+proj=longlat +ellps=WGS84 +no_defs")) prism_data <- data.frame() for (i in 1:12) { ppt_pd <- prism_archive_subset("ppt", "monthly normals", resolution = "800m", mon = i) ppt_pd <- pd_to_file(ppt_pd) ppt_pd_rast <- raster(ppt_pd) ppt <- raster::extract(ppt_pd_rast, points.spdf, fun=mean, na.rm=TRUE, sp=FALSE) tmean_pd <- prism_archive_subset("tmean", "monthly normals", resolution = "800m", mon = i) tmean_pd <- pd_to_file(tmean_pd) tmean_pd_rast <- raster(tmean_pd) tmean <- raster::extract(tmean_pd_rast, points.spdf, fun=mean, na.rm=TRUE, sp=FALSE) df <- data.frame(Site = points$Site, Lat = points$Lat, Lon = points$Lon, mon = i, ppt = ppt, tmean = tmean) prism_data <- rbind(prism_data, df) }
5.5. Water Balance Calculations
result <- prism_data %>% left_join(points, by = c("Site" = "Site")) %>% mutate(Lat = Lat.x, Lon = Lon.x) %>% group_by(Site) %>% arrange(mon) %>% # defaults mutate(slope = SLOPE_QGIS, aspect = ASPECT_QGIS, soil_max = 100, hock = 4) %>% mutate(F = get_f(tmean), RAIN = get_rain(ppt, F), SNOW = get_snow(ppt, F), PACK = get_pack(ppt, F), MELT = get_melt(SNOW, PACK, F), W = RAIN + MELT, Days = days_in_month(mon), DL = get_dl(mon, Days, Lat), A = abs(180 - abs(aspect - 225)), # folded aspect HL = get_hl(Lat, slope, A), e = 0.611 * exp((17.3 * tmean) / (tmean + 237.3)), PET = 29.8 * Days * DL * HL * (e / (tmean + 273.2)), #PET = 29.8 * Days * DL * (e / (tmean + 273.2))) %>% SOIL = get_soil(soil_max, W, PET), dSOIL = get_d_soil(SOIL), AET = get_aet(PET, dSOIL, W), D = PET - AET)
5.6. Figures
5.6.1. Helper Functions
mt_list <- c("PANAMINT RANGE", "WAH WAH MOUNTAINS", "SILVER PEAK RANGE", "HOT CREEK RANGE", "SNAKE RANGE", "SNAKE RANGE(MORIAH)")
5.6.2. AET x CWD
result %>% group_by(Site) %>% filter(Elev_m != TRUE) %>% summarise(D = sum(D), AET = sum(AET), T = mean(tmean), ppt = sum(ppt), elev = mean(Elev_m)) %>% left_join(points, by = "Site") %>% ggplot(mapping = aes(x = D, y = AET, color = MT_range)) + geom_point() + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE) + ggtitle("AET x CWD Hamon PRISM Lutz")
5.6.3. AET, CWD, PPT Plot
Re-create figure 4 from (Bentz et al. 2022)
data_long <- result %>% group_by(Site) %>% summarise(D = sum(D), AET = sum(AET), T = mean(tmean), ppt = sum(ppt), elev = mean(Elev_m), MT_range = MT_range) %>% distinct() %>% pivot_longer(c("D", "AET", "ppt")) data_long %>% transform(name=factor(name,levels=c("D", "ppt", "AET"))) %>% ggplot() + geom_point(mapping = aes(x = T, y = value, color = MT_range)) + gghighlight(MT_range %in% mt_list, use_direct_label = FALSE, calculate_per_facet = TRUE) + facet_wrap(~ name, ncol = 1, scales = "free")
6. Hypotheses
- Low AET values -> carbon starvation -> increased susceptibility to MPB
- High treeline temperature growth release. High temps -> accelerated growth rates in T-limited systems -> decreased wood density (?) -> increased susceptibility to MPB
- High temperature release reported in PILO in White Mountains by (Salzer et al. 2009)
- Beetle Hypothesis: Exploding beetle populations, enhanced by climate change, are overwhelming usual food source, PIFL. The beetles are forced to seek out novel food sources (i.e., PILO, which has previously been shown to be resistant to MPB (Bentz et al. 2017)). Beetle feeding could occur at low, undetected levels in PILO and not cause apparent harm, might be now passing threshold that can kill trees.
- PIFL promixity theory: This is the theory proposed by (Bentz et al. 2022). It is proposed that the presence of nearby infested PIFL led to the infestation of normally resistant PILO trees.
7. Questions
All of the models are "relative" and don't report absolute magnitudes that reflect conditions on the ground for all of the outputs. Snow might be closest, but soil moisture, evapotranspiration are likely to be different from reality. This is partly because PET regardless of method is for a "reference crop" or short sward of grass like a golf course rough. For this reason, we really need response data to test against your different results to understand which model is giving the most useful results. For example, even tho Penman PET is the gold standard physical basis for ET we found Hamon performs better when we use it in a water balance model to estimate spring flow in the desert and Oudin is better correlated with NDVI than Penman, also in desert settings. It looks like the simple water balance is generating results that are at least consistent with the unusual dryness of the Panamint range. Is that right? You outlined a few hypotheses. Which do you think is most plausible? What more could be done to evaluate each? What additional data would you need to resolve which hypothesis seems most reasonable? After making it this far, what do you think is the best course of action if were able to pick this project up again?
- The simple water balance results indicate that the Panamint range is exceptionally dry in certain measures, PPT, SOIL, and PACK. However, it is middle of the road among the points in the Bentz data set in CWD (As Bentz et al. found). In fact, all 4 of the ranges with reported PILO mortality had similar results for these measures in the simple water balance results. I'm not totally sure how to interpret these results, as my understanding is that CWD should indicate the dryness that the plants are experiencing, however if a measure like snowpack is so low, I would imagine it would have some physiological consequences.
8. Next Steps
8.1. TODO Validate R implementation of Simple Water Balance Model
Need to get calculations to match up with excel version, calculations are currently slightly off, although patterns are similar. T and P values are slightly off so there may be an issue with my daymet data aggregation code.
8.2. Field Work
- Gather information on sympatric PIFL trees. What are mortality and beetle infestation rates in PIFL, and how do they compare to PILO? Is there a threshold we can identify where PIFL experiences MPB infestation/mortality but PILO does not?
- Determine cause of mortality in sites not attributed to MPB - Silver Peak and Hot Creek Range
9. Discussion
9.1. Issues
Penman PET calculation returns extreme values for CWD and AET. This occurs in both analyses above implementing the Simple Water Balance model in R and Excel. The CWD values generated would be expected from a desert rather than montane forest(Tercek et al. 2021). This occurs with both our version as well as the version from WaterBalance. The patterns remain consistent between Hamon and Penman analyses, but absolute values vary. The differences in calculated CWD appear to result from the soil moisture content calculations (SOIL), with Penman calculations showing reduced availability of late season soil moisture compared with Hamon calculations.
Possible Reasons for differences between models:
- Additional Variables accounted for in Simple Water Balance Model that are missing from Lutz 2010 method
- jtemp - Lutz method uses simpler method to determine F
- Use of Hamon/Penman calculations
- vp - not used in
WaterBalance:ET_Hamon_daily - srad - not used in
WaterBalance:ET_Hamon_daily - SOIL/PACK - Penman/Hamon calculations produce significantly different results for SOIL/PACK. Within the Simple Water Balance model comparisons, this can account for the difference in magnitudes in CWD calculations, but this does not account for the differences in patterns betwen Hamon SWB and Hamon Lutz calculations.
- vp - not used in
- Calculations from daily time series (SWB) instead of 30 yr monthly normals (Lutz)
9.2. Conclusions
(Bentz et al. 2022) determined that the PILO mortality event on Telescope Peak occured at a location that is middle of the road in PILO's climate space. Their climatic analysis used PRISM 30 year normals and the Hamon method to determine PET, following (Lutz, van Wagtendonk, and Franklin 2010). While their code and thus exact method was not provided, we were able to approximately recreate their analysis above, receiving similar results.
Our analysis using the Simple Water Balance model method placed the Telescope Peak sites in a different climatic space than reported by (Bentz et al. 2022). Results varied based on choice of PET calculation method (Hamon vs Penman), with Penman returning much higher CWD and AET values than Hamon, and thus a hotter and drier environment. However, both PET calculations produce the same patterns and relative placement of locations within the Bristlecone climatic space. With both Penman and Hamon PET calculations, the Simple Water Balance model places the Telescope Peak points at the low extreme of AET values for all sites, and around average for CWD.
Additionally, the Simple Water Balance model shows low snow pack (PACK) values in the years preceding 2017 for the Telescope Peak points. Interestingly, (Hankin and Bisbing 2021) found improved PILO regeneration with decreasing spring snowpack, due to release from energy limitations and longer growing seasons. While possibly beneficial to PILO seedlings, this temperature release, also reported by (Salzer et al. 2009), may have deleterious effects on PILO or at least on their ability to resist MPB infestation.