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) )

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)
# --- File ---
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 (wait until pressure data or use - THIS IS RELATIVE TO OUR ZERO HEAD 

# --- 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
r = 0.01 #Radius of the well casing in m  
F_Shape = 2*pi*L/(log(L/D+ sqrt(1+(L/D)^2)))
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)
reading_start <- 1
reading_end   <- 160

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

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

dat <- dat %>% mutate(h_m = h_cm / 100)

# H0 = water level at t = 0 (first reading of the selected segment)
H0_m <- first(dat$h_m)

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   <- should be ~166 cm\n\n")
## Driving head H0 - H: 166.1297 cm   <- should be ~166 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
# 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()

The values make sense - The F shape is a little dubious we need to defend that better - and I used AI for the coding so I need to go through it in more detail but I think it mainly is okay. I will go through in more detail, but I thihk it is within the realm of reason.