This document is an analysis of karst water geochemistry collected in Cornwell Cave, Preston County, WV with an Onset Hobo U24-001 Conductivity Logger. Data was collected over a three-month period during summer of 2025, from May 18 - August 16.
Water infiltrates the ceiling of Cornwell’s 50-foot-tall “Rain Room”, having breached the overlying clastics at the Savage Dam and Loyalhanna contact, percolates through autochthonous clastic breakdown sediments on the floor, and then collects as a shallow stream which flows down-dip at the base of the Greenbrier Group. The data logger was placed just downstream of the “Rain Room” at the fully-formed stream, before “The Canyon”.
The source of the infiltrating water is unknown, but may be pirated from Hackelbarney Run which crosses the steep, narrow stripe outcropping of the Greenbrier Group 0.75 miles east of the cave, or may be originating from groundwater percolating from the walls of the Cheat River Canyon above the cave.
Data Import and Preparation
Read the CSV data as offloaded from Hobo U24 and exported from Hoboware 3.7.28. Discard unnecessary columns.
data <- read_csv("Cornwell_conductivity_2025.csv",
skip = 1, comment = "#",
col_names = c("n", "DateTime", "Conductivity_uS_cm", "Temperature_F")) |>
select(DateTime, Conductivity_uS_cm, Temperature_F) |>
mutate(DateTime = mdy_hms(DateTime)) # TODO: TZ should be EDT, this interprets them as UTC
summary(data)## DateTime Conductivity_uS_cm Temperature_F
## Min. :2025-05-18 00:00:00 Min. : 12.30 Min. :48.78
## 1st Qu.:2025-06-09 16:03:45 1st Qu.: 84.20 1st Qu.:49.24
## Median :2025-07-02 08:07:30 Median : 91.80 Median :49.55
## Mean :2025-07-02 08:07:30 Mean : 89.05 Mean :49.76
## 3rd Qu.:2025-07-25 00:11:15 3rd Qu.:105.50 3rd Qu.:49.98
## Max. :2025-08-16 16:15:00 Max. :139.30 Max. :53.17Calculate Specific Conductance
The U24 records raw conductivity values. Convert to Specific Conductance (K25) to compensate for temperature variance.
The term “specific conductance” is defined as the electrical conductance of 1 cm³ of a solution at 25 °C. If the electrical conductance is measured at another temperature, the value is corrected to what it would be at 25 °C and reported as specific conductance at 25 °C.
See USGS Water-Quality Data Manual: https://pubs.usgs.gov/tm/09/a6.3/tm9-a6_3.pdf
# Specific Conductance "K25" is the conductivity normalized to what it'd be at 25°C
# See USGS Water-Quality Data Manual: https://pubs.usgs.gov/tm/09/a6.3/tm9-a6_3.pdf
#
# I'm using simplified linear coefficient of 0.020
specific_conductance <- function(Conductivity_uS_cm, Temperature_C, TempCoefficient = 0.020) {
Conductivity_uS_cm / (1 + TempCoefficient * (Temperature_C - 25.0))
}
data <- data |>
mutate(
Temperature_C = (Temperature_F - 32) / 1.8,
SpecificConductance_uS_cm = specific_conductance(Conductivity_uS_cm, Temperature_C)
)
head(data)## # A tibble: 6 × 5
## DateTime Conductivity_uS_cm Temperature_F Temperature_C
## <dttm> <dbl> <dbl> <dbl>
## 1 2025-05-18 00:00:00 85.5 48.9 9.4
## 2 2025-05-18 00:15:00 85.5 48.9 9.4
## 3 2025-05-18 00:30:00 85.5 48.9 9.4
## 4 2025-05-18 00:45:00 85.5 48.9 9.4
## 5 2025-05-18 01:00:00 85.5 48.9 9.4
## 6 2025-05-18 01:15:00 85.5 48.9 9.4
## # ℹ 1 more variable: SpecificConductance_uS_cm <dbl>Join Recharge Events
The start of presumed precipitation events were manually
identified by patterns in the temperature and conductivity data and
saved as CSV. Join them here by adding an Event column to
associate each measurement with its precip event.
NOTE: No actual precipitation data has been correlated with this data yet!
# read manually-identified precip event datetimes
events <- read_csv("events.csv")
# add a column to represent the precip event index
data <- data |>
mutate(Event = cut(DateTime,
breaks = c(events$DateTime, max(data$DateTime) + days(1)),
labels = events$N,
right = FALSE,
include.lowest = TRUE))
# remove the first event 0 since it's an artifact of our logger deployment (we don't have its start)
events <- events |> filter(N != 0)
head(data)## # A tibble: 6 × 6
## DateTime Conductivity_uS_cm Temperature_F Temperature_C
## <dttm> <dbl> <dbl> <dbl>
## 1 2025-05-18 00:00:00 85.5 48.9 9.4
## 2 2025-05-18 00:15:00 85.5 48.9 9.4
## 3 2025-05-18 00:30:00 85.5 48.9 9.4
## 4 2025-05-18 00:45:00 85.5 48.9 9.4
## 5 2025-05-18 01:00:00 85.5 48.9 9.4
## 6 2025-05-18 01:15:00 85.5 48.9 9.4
## # ℹ 2 more variables: SpecificConductance_uS_cm <dbl>, Event <fct>Calculate Ratio
Join in the ratio of temperature to conductance to try and identify a baseline linear and anomalous nonlinear relationship.
# calculate the ratio of temp to conductance to find potential anomalies
data <- data |>
mutate(Ratio = Temperature_F / SpecificConductance_uS_cm)Write Data
Write updated data to CSV for external analysis.
write_csv(data, "Cornwell_conductivity_2025_updated.csv")Data Analysis
Time Series
Plot conductivity and temperature together to visualize our water parameters.
# Get ranges for scaling
cond_range <- range(data$SpecificConductance_uS_cm, na.rm = TRUE)
temp_range <- range(data$Temperature_F, na.rm = TRUE)
# Create scaling factor to fit temperature on conductivity scale
scale_factor <- diff(cond_range) / diff(temp_range)
temp_scaled <- (data$Temperature_F - temp_range[1]) * scale_factor + cond_range[1]
# Create a dual-axis plot (Conductivity and Temperature sharing X axis with independent Y)
p <- ggplot(data, aes(x = DateTime)) +
geom_vline(data = events, aes(xintercept = DateTime), color = "red", alpha = 0.4, linetype = "dashed") +
geom_line(aes(y = SpecificConductance_uS_cm), color = "black", linewidth = 0.8) +
geom_line(aes(y = temp_scaled), color = "blue", linewidth = 0.8) +
# Primary y-axis (left) - Conductivity
scale_y_continuous(
name = "Specific Conductance K25 (μS/cm)",
# Secondary y-axis (right) - Temperature
sec.axis = sec_axis(
trans = ~ (. - cond_range[1]) / scale_factor + temp_range[1],
name = "Temperature (°F)",
labels = function(x) paste0(round(x, 1), "°F")
)
) +
# X-axis - DateTime
scale_x_datetime(
name = NULL,
date_labels = "%Y-%m-%d",
date_breaks = "1 week"
) +
theme_minimal() +
theme(
axis.title.y.left = element_text(color = "black"),
axis.text.y.left = element_text(color = "black"),
axis.title.y.right = element_text(color = "blue"),
axis.text.y.right = element_text(color = "blue"),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5)
) +
labs(
title = "Cornwell Cave: Water Conductivity and Temperature"
)
print(p)
Ratio Analysis
Analyze ratio of temperature to conductivity to try and find anomalies
# Plot 3a: Simple ratios over time
p3a <- ggplot(data, aes(x = DateTime)) +
geom_vline(data = events, aes(xintercept = DateTime), color = "red", alpha = 0.4, linetype = "dashed") +
geom_line(aes(y = Ratio), linewidth = 0.8) +
scale_x_datetime(
name = NULL,
date_labels = "%Y-%m-%d",
date_breaks = "1 week"
) +
labs(
title = "Cornwell Cave: Temperature-to-Conductivity Ratio",
subtitle = "Higher values = temperature high relative to conductivity",
y = "Ratio (°F per μS/cm)",
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
print(p3a)
Conductivity-Temperature Relationship
Analyze the relationship between conductivity and temperature to show their non-linear relationship over the entire 3-month period. We use color to try and illustrate change over time.
# Analyze relationship between conductivity and temperature
# Scatter plot: Temperature vs Specific Conductance (colored by time)
p2 <- ggplot(data, aes(x = Temperature_F, y = SpecificConductance_uS_cm)) +
geom_point(aes(color = DateTime), alpha = 0.7, size = 1.2) +
scale_color_viridis_c(name = NULL, labels = function(x) format(as.POSIXct(x, origin = "1970-01-01"), "%Y-%m-%d")) +
labs(
title = "Specific Conductance (K25) vs Temperature",
x = "Temperature (°F)",
y = "Specific Conductance K25 (μS/cm)"
) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
print(p2)
Hysteresis Loop Analysis
For every precipitation event, analyze the non-linear hysteresis pattern as conductivity changes given an influx of precipitation into the karst system. Temperature is serving as a semi-proxy for discharge since we don’t have this data available. The assumption is that increasing temperature implies meteoric recharge from precipitation; the magnitude of this temperature delta is greater in the hotter late-summer months than the earlier start of our deployment.
library(patchwork)
# Get global scales for consistent axes across all plots
global_cond_range <- range(data$SpecificConductance_uS_cm, na.rm = TRUE)
global_temp_range <- range(data$Temperature_F, na.rm = TRUE)
global_time_range <- range(data$DateTime, na.rm = TRUE)
scale_factor <- diff(global_cond_range) / diff(global_temp_range)
# Create plots for each precip event
event_plots <- map(events$N, function(event_num) {
event_data <- data |> filter(Event == event_num)
event_datetime <- events$DateTime[events$N == event_num]
event_title <- paste0("Event ", event_num, ": ", format(event_datetime, "%Y-%m-%d %H:%M"))
# Time series plot with dual axis
p_timeseries <- ggplot(event_data, aes(x = DateTime)) +
geom_vline(xintercept = event_datetime, color = "red", linetype = "dashed", linewidth = 1) +
geom_line(aes(y = SpecificConductance_uS_cm), color = "black", linewidth = 0.8) +
geom_line(aes(y = (Temperature_F - global_temp_range[1]) * scale_factor + global_cond_range[1]),
color = "blue", linewidth = 0.8) +
scale_x_datetime(limits = global_time_range) +
scale_y_continuous(
name = "K25 (μS/cm)",
limits = global_cond_range,
sec.axis = sec_axis(
trans = ~ (. - global_cond_range[1]) / scale_factor + global_temp_range[1],
name = "Temperature (°F)"
)
) +
ggtitle(event_title) +
theme_minimal() +
theme(
axis.title.y.left = element_text(color = "black"),
axis.text.y.left = element_text(color = "black"),
axis.title.y.right = element_text(color = "blue"),
axis.text.y.right = element_text(color = "blue"),
plot.title = element_text(size = 10, hjust = 0.5, face = "bold")
) +
labs(x = NULL)
# Hysteresis loop plot
p_hysteresis <- ggplot(event_data, aes(x = Temperature_F, y = SpecificConductance_uS_cm)) +
geom_point(aes(color = DateTime), alpha = 0.7, size = 0.8) +
geom_point(data = slice_head(event_data, n = 1), shape = 17, size = 1.5, color = "red") +
scale_color_viridis_c(name = NULL, labels = function(x) format(as.POSIXct(x, origin = "1970-01-01"), "%m-%d")) +
scale_x_continuous(limits = global_temp_range) +
scale_y_continuous(limits = global_cond_range) +
ggtitle(paste("Hysteresis Loop", event_num)) +
labs(x = "Temperature (°F)", y = "K25 (μS/cm)") +
theme_minimal() +
theme(
plot.title = element_text(size = 10, hjust = 0.5, face = "bold"),
axis.title.x = element_text(color = "blue"),
axis.text.x = element_text(color = "blue")
)
# Combine side by side
p_timeseries + p_hysteresis
})
# Stack all event plot pairs vertically
reduce(event_plots, `/`) +
plot_annotation(
title = "Cornwell Cave Recharge Event Hysteresis Analysis",
theme = theme(plot.title = element_text(hjust = 0.5))
)