How dtc relates to vpd and swc
Derek Kober
2026-04-08
Brief
How abiotic factors relate to dTc. - SWC - ET - VPD
Created
2026-04-08 10:54:51.380201
Last Rendered
2026-04-08
TODO
Set Up
Load Libraries
Set Global Variables
Load Data
Soil
soil |>
ggplot(aes(x = ts, y = swc_mean, color = depth_label)) +
geom_line(alpha = 0.7, linewidth = 0.4) +
labs(
x = "Time",
y = "SWC (%)",
color = "Depth (cm)",
title = "SWC averaged across sensors at each depth (covered pits)"
) +
theme(
strip.text = element_text(size = 12),
legend.position = "bottom"
)
p_par <- ggplot(flux, aes(x = ts, y = par)) +
geom_line()
ggplotly(p_par)p_dtc_30 <- plot_ly(dtc_long) |>
add_lines(
x = ~ts,
y = ~dtc,
color = ~roi,
name = ~roi
) |>
layout(
yaxis = list(title = "dTc (°C)")
)
p_swc_30 <- plot_ly(soil) |>
add_lines(
x = ~ts,
y = ~swc_mean,
color = ~depth_label,
name = ~depth_label
) |>
layout(
yaxis = list(title = "SWC (%)")
)
p_vpd_30 <- plot_ly(flux) |>
add_lines(
x = ~ts,
y = ~vpd,
name = "VPD (30 min)"
) |>
layout(yaxis = list(title = "VPD (kPa)"))
p_et_30 <- plot_ly(flux) |>
add_lines(
x = ~ts,
y = ~et,
name = "ET (30 min)"
) |>
layout(yaxis = list(title = "ET (kg m-2 s-1)"))
p_par_30 <- plot_ly(flux) |>
add_lines(
x = ~ts,
y = ~par,
name = "PAR (30 min)"
) |>
layout(yaxis = list(title = "PAR (µmol m-2 s-1)"))
p_precip_30 <- plot_ly(flux) |>
add_lines(
x = ~ts,
y = ~precip,
name = "Precip"
) |>
layout(yaxis = list(title = "Precip (mm)"))
# CHANGED: include precip in the subplot because the title and notes both frame
# rain pulses as part of the interpretation.
fig_30min <- subplot(
p_dtc_30,
p_et_30,
p_par_30,
p_vpd_30,
p_swc_30,
p_precip_30,
nrows = 6,
shareX = TRUE,
titleY = TRUE
) |>
layout(
xaxis = list(title = "Time"),
title = "30-min dTc, SWC, VPD, ET, and Precip"
)## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colorsfig_30minOne important thing I noticed from the 30-min time series:
When the site is dry, ET tends to peak early in the morning (7-8am).
After rain events, ET peaks later, sometimes close to midday (9am-11am).
This shifting peak timing suggests we need to be careful about aligning ET with dTc, since dTc reliably peaks around midday. This shift might influence relationships under wet vs dry conditions.
Pre daily
In this section I’m trying to figure out how to aggregate climate variables—ET, VPD, and SWC—so they can be compared meaningfully with daily midday dTc. The main question is:
Should I compare dTc to full-day values, daytime values (PAR-based), or midday values (10:00–14:00)?
For ET, this choice matters because ET has strong diurnal patterns and the timing of the ET peak can shift depending on soil moisture. Currently i define daytime is defined using a PAR threshold >1. Midday is defined the same way as in the dTc calculations (10:00–14:00).
The goal here is to run a quick diagnostic and see which aggregation window (full day, daytime, or midday) is best aligned with daily dTc patterns.
I do the same for VPD, since VPD also has a pronounced diurnal cycle.
For both ET and VPD i set up time-series and pairwise plots to compare each aggregation against midday dTc.
# CHANGED: ET summaries are prepared in the load chunk now, so this section can
# focus on comparing aggregation windows instead of recalculating them.
p_et_total <- make_plotly_daily_metric(
daily_et,
y_var = "et_total_day",
series_name = "Total ET (24h)",
y_title = LAB_ET_DAILY
)
p_et_daytime <- make_plotly_daily_metric(
daily_et,
y_var = "et_total_daytime",
series_name = "Total ET (daytime)",
y_title = LAB_ET_DAILY
)
p_et_midday <- make_plotly_daily_metric(
daily_et,
y_var = "et_total_midday",
series_name = "Total ET (midday)",
y_title = LAB_ET_DAILY
)
fig_et <- subplot(
p_et_30,
p_et_total,
p_et_daytime,
p_et_midday,
nrows = 4,
shareX = TRUE,
titleY = TRUE
) |>
layout(
xaxis = list(title = "Date"),
title = "Comparison of ET daily totals"
)
fig_etet_dtc <- daily_et |>
left_join(dtc_midday, by = "date") |>
drop_na()
et_dtc_long <- et_dtc |>
pivot_longer(
cols = c(et_total_day, et_total_daytime, et_total_midday),
names_to = "et_agg",
values_to = "et"
)
et_dtc_aug <- et_dtc_long |>
augment_species_daily_fit(x_var = "et", agg_var = "et_agg")
p_et_species <- plot_agg_vs_dtc(
data = et_dtc_aug,
x_var = "et",
agg_var = "et_agg",
x_label = "Daily ET (kg m-2 day-1)",
title_text = "Midday dTc vs ET \nROIs colored; species-level regression in black"
)
p_et_species## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Notes on ET
Time series of total ET over 24 hours and total ET daytime look very similar. The daytime version is slightly lower in magnitude but follows almost the exact same pattern. Midday ET (10:00–14:00), on the other hand, is much coarser and naturally smaller in magnitude because it only covers a tight window, though the overall daily trend is still visible.
# CHANGED: use the same helper as ET so the VPD aggregation comparison stays in
# the same format and only needs one plotting style to maintain.
p_vpd_day <- make_plotly_daily_metric(
daily_vpd,
y_var = "vpd_mean_day",
series_name = "24h Mean VPD",
y_title = "VPD (kPa)"
)
p_vpd_daytime <- make_plotly_daily_metric(
daily_vpd,
y_var = "vpd_mean_daytime",
series_name = "Daytime Mean VPD",
y_title = "VPD (kPa)"
)
p_vpd_midday <- make_plotly_daily_metric(
daily_vpd,
y_var = "vpd_mean_midday",
series_name = "Midday Mean VPD",
y_title = "VPD (kPa)"
)
fig_vpd <- subplot(
p_vpd_30,
p_vpd_day,
p_vpd_daytime,
p_vpd_midday,
nrows = 4,
shareX = TRUE,
titleY = TRUE
) |>
layout(
xaxis = list(title = "Date"),
title = "Comparison of VPD Means by Aggregation Window"
)
fig_vpd# Join VPD aggregations with midday dTc
vpd_dtc <- daily_vpd |>
left_join(dtc_midday, by = "date") |>
drop_na()
vpd_dtc_long <- vpd_dtc |>
pivot_longer(
cols = starts_with("vpd_mean"),
names_to = "vpd_agg",
values_to = "vpd"
)
vpd_dtc_aug <- vpd_dtc_long |>
augment_species_daily_fit(x_var = "vpd", agg_var = "vpd_agg")
p_vpd_species <- plot_agg_vs_dtc(
data = vpd_dtc_aug,
x_var = "vpd",
agg_var = "vpd_agg",
x_label = "VPD (kPa)",
title_text = "Midday dTc vs VPD Aggregations\nIndividual ROIs Colored; Species-Level Pooled Regression in Black"
)
p_vpd_species## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Notes on VPD
The midday mean VPD is slightly smoother with fewer fluctuations, but the overall patterns across the three are nearly the same.
Pre Daily summary
Notes on Choosing Daily Aggregation Windows for ET, VPD, SWC, and dTc Comparisons
Working through the daily comparisons of evapotranspiration (ET), VPD, soil water content (SWC), and midday dTc. I want to be consistent about which daily aggregation makes the most sense for each variable.
- ET: For ET, I’m leaning toward using total daytime ET rather than total 24-hour or midday-only. Right now, “daytime” is defined using a PAR threshold of 2. I haven’t fully checked that threshold yet, but visually it works: PAR stays below ~1 at night, spikes above 1 quickly after sunrise, and using 2 helps buffer noise. When looking at ET patterns at the flux tower: Under dry conditions, the ET peak tends to occur early in the morning (e.g., June 20: peak around 7–8 AM). After rain events, ET peaks shift later (e.g., June 25: peaks closer to ~9:30–11 AM). Meanwhile, dTc reliably peaks around midday (noon–1 PM). Because ET and dTc peaks often don’t line up in time, using the midday ET would miss early-morning ET peaks on dry days. So, to avoid misalignment issues, I think total daytime ET is the best metric to compare with midday dTc.
- VPD: For VPD, I’m less certain. VPD and dTc tend to peak at similar times of day, so daytime and midday averages both correlate well with dTc. In my pairwise comparisons, both are significant and explain more variance than a full-day average. I’ll check with Marcy whether daytime or midday makes more sense, but at first glance either seems fine.
- SWC: For SWC, I’ll continue using the full-day average. SWC doesn’t show the strong diurnal dynamics that ET or VPD do, so daily mean SWC is okay and avoids adding unnecessary complexity.
Daily
# CHANGED: reuse the daily tables created in the load chunk so this summary
# figure stays aligned with the pairwise comparisons below.
p_daily_dtc <- plot_ly(dtc_midday) |>
add_lines(
x = ~date,
y = ~dtc_midday,
color = ~roi,
name = ~roi,
line = list()
) |>
add_markers(
x = ~date,
y = ~dtc_midday,
color = ~roi, # <-- ensures matching colors
marker = list(size = 5),
showlegend = FALSE # hides markers from legend
) |>
layout(
yaxis = list(title = "Midday dTc (°C)")
)
p_daily_swc <- plot_ly(daily_swc) |>
add_lines(
x = ~date,
y = ~swc_day_mean,
color = ~depth_label,
name = ~depth_label
) |>
add_markers(
x = ~date,
y = ~swc_day_mean,
marker = list(size = 3, color = "black"),
showlegend = FALSE
) |>
layout(
yaxis = list(title = "SWC (%)")
)
p_precip_daily <- plot_ly(daily_precip) |>
add_bars(
x = ~date,
y = ~precip_total,
name = "Daily Precip",
marker = list(line = list(width = 0))
) |>
layout(
yaxis = list(title = "Total Daily Precip (mm)"),
xaxis = list(title = "Date"),
title = "Precipitation",
bargap = 0.1
)
fig_daily <- subplot(
p_daily_dtc,
p_daily_swc,
p_et_daytime,
p_vpd_daytime,
p_precip_daily,
nrows = 5,
shareX = TRUE,
titleY = TRUE
) |>
layout(
xaxis = list(title = "Date"),
title = "Daily Midday dTc vs Daily Mean SWC"
)## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(max(N, 3L), "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colorsfig_daily# CHANGED: keep the paired-response plotting logic in one helper and let SWC
# reuse it from a long table so depth can be faceted in one condensed plot.
pair_plot <- function(
df,
xvar,
xlab,
group_vars = character(),
facet_var = NULL,
facet_ncol = 1
) {
df_all <- df |>
filter(dtc_midday < 15) |>
bind_rows(
df |>
group_by(across(all_of(c(group_vars, "species", "date")))) |>
summarise(
x = mean(.data[[xvar]], na.rm = TRUE),
dtc_midday = mean(dtc_midday, na.rm = TRUE),
roi = "all",
.groups = "drop"
)
) |>
mutate(x = ifelse(roi == "all", x, .data[[xvar]]))
# --- NEW COLOR PALETTE LOGIC ---
# 1. Grab all the actual ROI names (excluding the 'all' group we just made)
actual_rois <- unique(df$roi)
# 2. Sort them into Piñon and Juniper arrays based on the first letter
p_names <- sort(grep("^p", actual_rois, ignore.case = TRUE, value = TRUE))
j_names <- sort(grep("^j", actual_rois, ignore.case = TRUE, value = TRUE))
# 3. Generate color gradients for however many of each tree you have
# Piñons get a gradient of light blue to dark blue
p_cols <- colorRampPalette(c("#9ecae1", "#084594"))(length(p_names))
# Junipers get a gradient of yellow to golden-orange (so they are visible!)
j_cols <- colorRampPalette(c("#fed976", "#cc4c02"))(length(j_names))
# 4. Tie the specific colors to the specific tree names
roi_pal <- setNames(c(p_cols, j_cols), c(p_names, j_names))
# -------------------------------
p <- ggplot() +
geom_point(
data = df_all |> filter(roi != "all"),
aes(x = x, y = dtc_midday, color = roi),
alpha = 0.5,
size = 0.7
) +
geom_smooth(
data = df_all |> filter(roi != "all"),
aes(x = x, y = dtc_midday, color = roi),
method = "lm", alpha = 0.5, se = FALSE, linewidth = 0.3
) +
geom_smooth(
data = df_all |> filter(roi == "all"),
aes(x = x, y = dtc_midday),
method = "lm", se = TRUE, color = "black", linewidth = 1.4
) +
stat_cor(
data = df_all |> filter(roi == "all"),
aes(x = x, y = dtc_midday),
method = "pearson",
label.x.npc = 0.97,
label.y.npc = 0.97,
hjust = 1,
size = 3
) +
# Add the generated palette to the plot!
scale_color_manual(values = roi_pal) +
labs(
x = xlab,
y = "Midday ∆T (°C)",
color = "ROI"
) +
theme_bw() +
theme(
strip.text = element_text(size = 10),
legend.position = "none"
) +
ylim(NA,15)
if (!is.null(facet_var)) {
p <- p +
facet_wrap(
vars(!!rlang::sym(facet_var)),
scales = "free_x",
ncol = facet_ncol,
strip.position = "right"
) +
theme(strip.text.y = element_text(angle = -90)) # Keeps side labels readable
}
p
}
swc_dtc_long <- swc_dtc |>
pivot_longer(
cols = starts_with("swc_"),
names_to = "depth_label",
names_pattern = "swc_(.*)_day_mean",
values_to = "swc_day_mean"
) |>
mutate(
depth_label = factor(depth_label, levels = c("5cm", "10cm", "30cm"))
)
p_swc_dtc_pairs <- pair_plot(
swc_dtc_long,
xvar = "swc_day_mean",
xlab = "Daily Mean SWC (%)",
group_vars = "depth_label",
facet_var = "depth_label"
)
p_vpd <- pair_plot(vpd_dtc, "vpd_mean_daytime", "Daytime Mean VPD (kPa)")
p_et <- pair_plot(et_dtc, "et_total_daytime", "Daytime Total ET (kg m-2 day-1)")
p_vpd_et_swc_dtc_pair <- ( p_vpd / p_et/p_swc_dtc_pairs) +
plot_annotation( title = "Pairwise SWC, VPD, and ET vs Midday dTc", tag_levels = "a") +
plot_layout(axis_titles = "collect")+
plot_layout(heights = c(1,1,4))
p_vpd_et_dtc_pair <- p_vpd/p_et +
plot_layout(axis_titles = "collect")
ggsave(filename = here("output/research_days/p_vpd_et_swc_dtc_pair.png"),p_vpd_et_swc_dtc_pair, width = 7, height = 10, units = "in",dpi = 300 )## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'#ggsave(filename = here("output/research_days/p_swc_dtc_pairs.png"),p_swc_dtc_pairs, width = 7, height = 5, units = "in" )Research day plots 2026
Time series of P and J smoothed? with precip and vpd and swc
plot_min_date <- min(dtc_midday$date, na.rm = TRUE)
plot_max_date <- max(dtc_midday$date, na.rm = TRUE)
dtc_midday <- dtc_midday |>
mutate(
species = factor(species, levels = c("piñon", "juniper"))
)
daily_swc <- daily_swc |>
filter(date >= plot_min_date, date <= plot_max_date) |>
mutate(
depth = factor(depth, levels = sort(unique(depth)))
)
# Main panel: midday ΔTc
p_dtc_midday <- dtc_midday |>
ggplot(aes(x = date, y = dtc_midday, color = species)) +
geom_point(alpha = 0.28, size = 1.2) +
geom_smooth(
method = "loess",
span = 0.25,
se = TRUE,
linewidth = 1.2
) +
labs(
y = expression("Midday ∆T (°C)"),
colour = "Species",
x = "Date"
) +
scale_x_date(
limits = c(plot_min_date, plot_max_date),
date_breaks = "1 month",
date_labels = "%b"
) +
theme_bw() +
scale_color_manual(values = SPECIES_COLS)
p_daily_et <- daily_et |> ggplot(aes(x = date, y = et_total_daytime)) + geom_line()
#----------------------------
# 2) Midday VPD context
#----------------------------
p_vpd <- daily_vpd |>
filter(date >= plot_min_date, date <= plot_max_date) |>
ggplot(aes(x = date, y = vpd_mean_midday)) +
geom_line(linewidth = 0.8) +
labs(
title = "Midday Vapor Pressure Deficit",
y = "VPD (kPa)",
x = "Date"
) +
scale_x_date(
limits = c(plot_min_date, plot_max_date),
date_breaks = "1 month",
date_labels = "%b"
) +
theme_bw()+
theme(legend.position = "none")
#----------------------------
# 3) Soil water content context
#----------------------------
max_precip <- max(daily_precip$precip_total, na.rm = TRUE)
max_swc <- max(daily_swc$swc_day_mean, na.rm = TRUE)
# If max_precip is 40 and max_swc is 20, factor is 2.
scale_factor <- max_precip / max_swc
# 2. Build the combined plot
p_combined_swc_precip <- ggplot() +
# Primary Axis (Left): Precipitation as columns
geom_col(
data = daily_precip |> filter(date >= plot_min_date, date <= plot_max_date),
aes(x = date, y = precip_total),
fill = "blue",
alpha = 0.7
) +
# Secondary Axis (Right): SWC as lines (scaled by the factor)
geom_line(
data = daily_swc,
aes(x = date, y = swc_day_mean * scale_factor, color = depth),
linewidth = 0.8
) +
# Define the Dual Y-Scales
scale_y_continuous(
name = "Precip. (mm)",
sec.axis = sec_axis(~ . / scale_factor, name = "SWC (%)")
) +
scale_x_date(
limits = c(plot_min_date, plot_max_date),
date_breaks = "1 month",
date_labels = "%b"
) +
# Using a nice color palette for soil depths (browns/earth tones often work well)
scale_color_viridis_d(option = "plasma", end = 0.8) +
labs(
x = "Date",
color = "Soil Depth (cm)"
) +
theme_bw()
p_combined_swc_precip## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_col()`).
#----------------------------
# Combine
#----------------------------
p_context <- (p_dtc_midday / p_vpd / p_combined_swc_precip) +
plot_layout(
guides = "collect",
axes = "collect"
) +
plot_annotation(tag_levels = 'a') &
labs(title = NULL, x = NULL) &
theme(legend.position = "bottom")
p_context## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_col()`).
ggplotly(p_dtc_midday)## `geom_smooth()` using formula = 'y ~ x'## Warning in is.na(u): is.na() applied to non-(list or vector) of type
## 'expression'#ggsave(filename = here("output/research_days/dtc_flux_context.png"), plot = p_context, height = 7, width = 10, units = "in", dpi = 300)# CHANGED: replace the old placeholder with a compact sanity check of the
# daily joined inputs used above so the chunk stays runnable and informative.
list(
swc_dtc = swc_dtc |>
select(date, species, roi, starts_with("swc_")) |>
head(),
vpd_dtc = vpd_dtc |>
select(date, species, roi, vpd_mean_daytime, dtc_midday) |>
head(),
et_dtc = et_dtc |>
select(date, species, roi, et_total_daytime, dtc_midday) |>
head()
)## $swc_dtc
## # A tibble: 6 × 6
## date species roi swc_5cm_day_mean swc_10cm_day_mean swc_30cm_day_mean
## <date> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2025-04-08 juniper j22 7.76 9.77 11.9
## 2 2025-04-08 juniper j23 7.76 9.77 11.9
## 3 2025-04-08 piñon p2 7.76 9.77 11.9
## 4 2025-04-08 piñon p21 7.76 9.77 11.9
## 5 2025-04-09 juniper j22 7.83 9.95 11.9
## 6 2025-04-09 juniper j23 7.83 9.95 11.9
##
## $vpd_dtc
## # A tibble: 6 × 5
## date species roi vpd_mean_daytime dtc_midday
## <date> <chr> <chr> <dbl> <dbl>
## 1 2025-04-08 juniper j22 1.48 5.24
## 2 2025-04-08 juniper j23 1.48 3.94
## 3 2025-04-08 piñon p2 1.48 4.08
## 4 2025-04-08 piñon p21 1.48 4.18
## 5 2025-04-09 juniper j22 1.85 1.22
## 6 2025-04-09 juniper j23 1.85 -0.523
##
## $et_dtc
## # A tibble: 6 × 5
## date species roi et_total_daytime dtc_midday
## <date> <chr> <chr> <dbl> <dbl>
## 1 2025-04-08 juniper j22 0.483 5.24
## 2 2025-04-08 juniper j23 0.483 3.94
## 3 2025-04-08 piñon p2 0.483 4.08
## 4 2025-04-08 piñon p21 0.483 4.18
## 5 2025-04-09 juniper j22 0.376 1.22
## 6 2025-04-09 juniper j23 0.376 -0.523# CHANGED: summarize the all/month/week diurnal curves through shared helpers so
# the SE calculations and plot styling stay consistent.
dtc_long <- dtc_long |>
mutate(
date = as_date(ts),
month = month(ts, label = TRUE, abbr = TRUE),
week = floor_date(ts, unit = "week", week_start = 1),
time_of_day = as_hms(format(ts, "%H:%M:%S"))
)
dtc_diurnal_all <- summarise_dtc_diurnal(dtc_long)
plot_dtc_diurnal(
data = dtc_diurnal_all,
title_text = "Mean diurnal ΔTc pattern by species",
breaks_vec = hms::as_hms(sprintf("%02d:00:00", seq(0, 23, by = 2)))
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
dtc_diurnal_month <- summarise_dtc_diurnal(dtc_long, group_vars = "month")
plot_dtc_diurnal(
data = dtc_diurnal_month,
title_text = "Mean diurnal ΔTc pattern by species and month",
facet_var = "month",
ncol = 3,
breaks_vec = hms::as_hms(sprintf("%02d:00:00", seq(0, 23, by = 4)))
)## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).
dtc_diurnal_week <- summarise_dtc_diurnal(dtc_long, group_vars = "week")
plot_dtc_diurnal(
data = dtc_diurnal_week,
title_text = "Mean diurnal ΔTc pattern by species and week",
facet_var = "week",
ncol = 2,
breaks_vec = hms::as_hms(c("00:00:00", "06:00:00", "12:00:00", "18:00:00"))
)## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_ribbon()`).