Ks is being measured using a direct- push Piezometer with falling head( adapted from Ronkanen & Klve, 2005) this program should be used with data from the device created by the FENCES project intended to measure KS in situ.

The Dimensions of our device:

Petronille Measured the Device:

The Dimensions are as follows:

Inner diameter of the water tank: 14.5 cm Tube Section (length: 109 cm, Ø 1.45 cm) Tube Section Ext 1 (length: 51 cm Ø = 1.20 cm) Tube Section Ext 2 (length: 49 cm Ø = 1.20cm) output is Ø = 1 cm and center of the circle is located at 5 cm from bottom

Theory - the rate of outflow (q at the peizometer tip at any (t) is porportional to the hydraulic conductivity of the soil and the unrecovered head difference H-h (Hvorselv 1951) )

What is unrecovered head difference ? -

H is the ground head h is the starting head (water level at top of reservoir) F = shape function K = ks - hydrological conductivity

q(t) = pi r ^2 dh/dt = FK(H-h)

That is the equation which we can re-organize to the following: as we do not know the q(t) directly but we do have the rate at which the reservoir decreases:

pi r ^2 dh/dt = FK(H-h)

re-arranging:

dh/dt + (FK/pir^2)h = (FK/pir^2)H

Solve:

ln((h-H)/(H0-H)) = -(FK/pir^2)*t

we can solve the right side of the equation - it is the slope of line

This was a bunch of theory but it is simple:

the slope of our line is the left side of the equation and then we can solve for k with a known F

# importing important libs

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.5.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
# importing file - note you have to edit the data file from the device - save it as a csv and get rid of the start headers 

csv_path <- "C:/Users/KAUFMADA/OneDrive - Trinity College Dublin/Data/Ks_measurements/Test_Ks_Measurement_Data_12_06_2026_R.csv"

datetime_col <- "datetime"         # name of datetime column in CSV
waterlevel_col <- "height_cm" # name of water level column in CSV (in cm)
# --- Device dimensions ---
tank_diameter_cm <- 14.5          # inner diameter of reservoir tank (cm)

# --- Field measurements ---
H_cm <- (-156)     # equilibrium ground water head before test (cm) 
                   # use the closest ph tube data
# ---------------------------------------------------------------------------------
# note we need convention here as the measurements depends on how we label things. For example if we say 0 is the zero in the reservoir then the ground water level is -x but starting water level is positive 20 if it is at the 20cm mark. The difference will always be the same its just what is easiest to type ? 
# ---------------------------------------------------------------------------------


# --- Shape factor ---
# F accounts for geometry of flow at piezometer tip
# From Hvorslev (1951) for a piezometer with L/R ratio
# You can enter a value from literature or calculate below
L = 0.01 #length of perforation on the tube in m
D = 0.0140 #diameter of the tube where there is the perforations in m
F_Shape = 2*pi*L/(log(L/D+ sqrt(1+(L/D)^2)))

# THE ONE THING HERE IS I MIGHT NEED TO UPDATE THIS AS the aspect ratio is not that large so we cannot simplify in the way I did ????????

cat("Calculated F from tip geometry:", round(F_Shape, 3), "\n")
## Calculated F from tip geometry: 0.095
dat <- read_csv(csv_path, show_col_types = FALSE)

# Falling-head segment only: readings 9-179
# (180+ is the disturbance/crash, not falling head)
# Before running program - look at the index in the csv of the measurement you want to focus on. 

reading_start <- 0
reading_end   <- 160

# creating new columns - copies of time and height for manipulation - the time we subtract to start the index at 0 

dat <- dat %>%
  filter(reading_n >= reading_start,
         reading_n <= reading_end) %>%
  mutate(h_cm  = height_cm,
         t_sec = reading_n - first(reading_n))   # each reading = 1 s
# Convert cm to metres - as Ks is measured in meters - simply dividing. 

r_m <- (tank_diameter_cm / 2) / 100    # reservoir radius (m): 14.5 cm tank -> 0.0725
H_m <- H_cm / 100                      # equilibrium ground head (m), SAME DATUM as h
                                       # the tank's zero (residual) level. AS MENTIONED ABOVE WE CAN CHOOSE WHAT WE WANT TO USE AS ZERO.

dat <- dat %>% mutate(h_m = h_cm / 100) # getting current head in meters.

# H0 = water level at t = 0 (first reading of the selected segment)
H0_m <- first(dat$h_m) # as the device measures from 0 in the reservoir this is cm above this 0 point. So + not - head. 

cat("H0 (initial reservoir head):", H0_m * 100, "cm\n")
## H0 (initial reservoir head): 10.1297 cm
cat("H  (equilibrium ground head):", H_cm, "cm\n")
## H  (equilibrium ground head): -156 cm
cat("Driving head H0 - H:", (H0_m - H_m) * 100, "cm")
## Driving head H0 - H: 166.1297 cm
# ln((h - H) / (H0 - H))
# With H = -156, both numerator and denominator are large positive numbers
# (~160 cm), so the ratio starts at 1 and shrinks slightly -> ln from 0 downward
dat <- dat %>%
  mutate(ln_term = log((h_m - H_m) / (H0_m - H_m)))

# Drop any non-finite values (can't occur here unless h reaches H, but keep as guard)
dat <- dat %>% filter(is.finite(ln_term))
# ============================================================
# SECTION 5: LINEAR REGRESSION TO GET SLOPE
# ============================================================
# Fit y = m*t  (no intercept - forced through origin because
# ln((h-H)/(H0-H)) = 0 at t = 0 by construction)

fit <- lm(ln_term ~ 0 + t_sec, data = dat)

slope_m   <- coef(fit)["t_sec"]   # this is -(F*K / (pi*r^2)), units 1/s
r_squared <- summary(fit)$r.squared

cat("Slope (m):", slope_m, "   <- expect ~ -1.2e-04 with H = -156\n")
## Slope (m): -0.000119642    <- expect ~ -1.2e-04 with H = -156
cat("R-squared:", round(r_squared, 4), "\n\n")
## R-squared: 0.9896
# plugging in our values

# K = -(slope * pi * r^2) / F
# slope in 1/s, r in m  ->  K in m/s

K_ms <- -(slope_m * pi * r_m^2) / F_Shape

cat("=============================================\n")
## =============================================
cat("Hydraulic Conductivity K =", formatC(K_ms, format = "e", digits = 3), "m/s\n")
## Hydraulic Conductivity K = 2.089e-05 m/s
cat("=============================================\n")
## =============================================
dat <- dat %>%
  mutate(ln_predicted = slope_m * t_sec)

ggplot(dat, aes(x = t_sec)) +
  geom_point(aes(y = ln_term), colour = "steelblue", size = 2) +
  geom_line(aes(y = ln_predicted), colour = "red", linewidth = 1) +
  labs(
    title = paste0("Falling Head Test - K = ",
                   formatC(K_ms, format = "e", digits = 3), " m/s"),
    subtitle = paste0("R\u00b2 = ", round(r_squared, 4),
                      "  |  F = ", round(F_Shape, 4), " m",
                      "  |  H = ", H_cm, " cm",
                      "  |  n = ", nrow(dat), " points"),
    x = "Elapsed time (seconds)",
    y = "ln((h - H) / (H0 - H))",
    caption = "Points = observed | Red line = no-intercept fit"
  ) +
  theme_minimal()

Lets do a little bit of investigating - does this behave how we expect it to ?

dat <- dat %>%
  mutate(
    h_fit_m   = H_m + (H0_m - H_m) * exp(slope_m * t_sec),
    h_fit_cm  = h_fit_m * 100,
    residual  = h_cm - h_fit_cm          # observed minus fitted (cm)
  )


ggplot(dat, aes(x = t_sec, y = h_cm)) +
  geom_line(colour = "steelblue", linewidth = 0.7) +
  geom_point(size = 1.2, colour = "steelblue") +
  labs(
    title = "Raw water level",
    x = "Elapsed time (s)", y = "h (cm)"
  ) +
  theme_minimal()

ggplot(dat, aes(x = t_sec, y = residual)) +
  geom_hline(yintercept = 0, colour = "red", linewidth = 0.8) +
  geom_point(colour = "steelblue", size = 1.8, alpha = 0.7) +
  geom_smooth(method = "loess", se = FALSE,
              colour = "orange", linewidth = 0.9, linetype = "dashed") +
  labs(
    title = "Residuals (observed \u2212 fitted) – should be random noise",
    x     = "Elapsed time (s)",
    y     = "Residual (cm)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'