diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index cac19cc..a0502d1 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -49,5 +49,6 @@ jobs:
           upload-snapshots: true
 
       - name: Test coverage
+        if: matrix.config.os == 'ubuntu-latest' && matrix.config.r == 'release'
         run: covr::codecov(quiet = FALSE)
         shell: Rscript {0}
diff --git a/tests/testthat/test_helper_functions.R b/tests/testthat/test_helper_functions.R
new file mode 100644
index 0000000..6185d8f
--- /dev/null
+++ b/tests/testthat/test_helper_functions.R
@@ -0,0 +1,135 @@
+## ---- has.vars ----
+
+test_that("has.vars returns TRUE for existing columns", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20, wnd = 5)
+  expect_true(has.vars(df, "wtr"))
+  expect_true(has.vars(df, "wnd"))
+})
+
+test_that("has.vars returns FALSE for missing columns", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20)
+  expect_false(has.vars(df, "wnd"))
+})
+
+test_that("has.vars is case-insensitive", {
+  df <- data.frame(datetime = Sys.time(), WTR = 20)
+  expect_true(has.vars(df, "wtr"))
+})
+
+test_that("has.vars handles multiple variable names", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20, wnd = 5)
+  result <- has.vars(df, c("wtr", "wnd", "par"))
+  expect_equal(result, c(TRUE, TRUE, FALSE))
+})
+
+test_that("has.vars errors when input is not a data.frame", {
+  expect_error(has.vars(c(1, 2, 3), "x"))
+})
+
+## ---- get.vars ----
+
+test_that("get.vars returns data.frame with datetime and matched columns", {
+  df <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01", "2020-01-02"), tz = "UTC"),
+    wtr = c(15, 16),
+    wnd = c(3, 4))
+  result <- get.vars(df, "wtr")
+  expect_s3_class(result, "data.frame")
+  expect_true("datetime" %in% names(result))
+  expect_true("wtr" %in% names(result))
+  expect_false("wnd" %in% names(result))
+})
+
+test_that("get.vars errors when no variable pattern matches", {
+  df <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wtr = 20)
+  expect_error(get.vars(df, "nonexistent_var"))
+})
+
+test_that("get.vars errors when no datetime column is present", {
+  df <- data.frame(wtr = 20, wnd = 5)
+  expect_error(get.vars(df, "wtr"))
+})
+
+## ---- rmv.vars ----
+
+test_that("rmv.vars removes named column", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20, wnd = 5)
+  result <- rmv.vars(df, "wnd")
+  expect_false("wnd" %in% names(result))
+  expect_true("wtr" %in% names(result))
+})
+
+test_that("rmv.vars with ignore.missing=TRUE does not error on absent column", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20)
+  expect_null(rmv.vars(df, "nonexistent", ignore.missing = TRUE))
+})
+
+test_that("rmv.vars with ignore.missing=FALSE errors on absent column", {
+  df <- data.frame(datetime = Sys.time(), wtr = 20)
+  expect_error(rmv.vars(df, "nonexistent", ignore.missing = FALSE))
+})
+
+## ---- get.offsets (from rLakeAnalyzer) ----
+
+test_that("get.offsets parses numeric suffix from column name", {
+  df <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wtr_1.5 = 20)
+  wtr_df <- get.vars(df, "wtr")
+  offsets <- get.offsets(wtr_df)
+  expect_equal(offsets, 1.5)
+})
+
+test_that("get.offsets returns NA for column without numeric suffix", {
+  df <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wtr = 20)
+  wtr_df <- get.vars(df, "wtr")
+  offsets <- suppressWarnings(get.offsets(wtr_df))
+  expect_true(is.na(offsets))
+})
+
+## ---- date2doy ----
+
+test_that("date2doy returns 1.5 for Jan 1 at noon UTC", {
+  dt <- as.POSIXct("2020-01-01 12:00:00", tz = "UTC")
+  result <- LakeMetabolizer:::date2doy(dt)
+  expect_equal(result, 1.5)
+})
+
+test_that("date2doy returns 1.0 for Jan 1 at midnight UTC", {
+  dt <- as.POSIXct("2020-01-01 00:00:00", tz = "UTC")
+  result <- LakeMetabolizer:::date2doy(dt)
+  expect_equal(result, 1.0)
+})
+
+test_that("date2doy handles vector input", {
+  dts <- as.POSIXct(c("2019-01-01 00:00:00", "2019-01-01 12:00:00"), tz = "UTC")
+  result <- LakeMetabolizer:::date2doy(dts)
+  expect_length(result, 2)
+  expect_equal(result, c(1.0, 1.5))
+})
+
+test_that("date2doy Dec 31 returns ~365 for non-leap year", {
+  dt <- as.POSIXct("2019-12-31 00:00:00", tz = "UTC")
+  result <- LakeMetabolizer:::date2doy(dt)
+  expect_equal(result, 365)
+})
+
+## ---- calc.freq ----
+
+test_that("calc.freq returns 144 for 10-minute interval data", {
+  start <- as.POSIXct("2020-01-01 00:00:00", tz = "UTC")
+  dt <- start + seq(0, 143 * 600, by = 600)  # 144 obs, 10-min apart
+  result <- LakeMetabolizer:::calc.freq(dt)
+  expect_equal(result, 144)
+})
+
+test_that("calc.freq returns 48 for 30-minute interval data", {
+  start <- as.POSIXct("2020-01-01 00:00:00", tz = "UTC")
+  dt <- start + seq(0, 47 * 1800, by = 1800)  # 48 obs, 30-min apart
+  result <- LakeMetabolizer:::calc.freq(dt)
+  expect_equal(result, 48)
+})
diff --git a/tests/testthat/test_k600.2.kGAS.R b/tests/testthat/test_k600.2.kGAS.R
new file mode 100644
index 0000000..663cbfa
--- /dev/null
+++ b/tests/testthat/test_k600.2.kGAS.R
@@ -0,0 +1,55 @@
+test_that("k600.2.kGAS.base formula: kGAS = k600 * (Sc/600)^-0.5", {
+  k600 <- 2.0
+  temp  <- 20
+  Sc    <- getSchmidt(temp, "O2")
+  expected <- k600 * (Sc / 600)^(-0.5)
+  result   <- k600.2.kGAS.base(k600, temp, "O2")
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("k600.2.kGAS.base: kGAS > k600 at warm water where Sc(O2) < 600", {
+  # At 20°C, Sc(O2) ≈ 531 < 600, so (Sc/600)^-0.5 > 1
+  k600 <- 2.0
+  kGAS <- k600.2.kGAS.base(k600, 20, "O2")
+  expect_true(kGAS > k600)
+})
+
+test_that("k600.2.kGAS.base: kGAS ≈ k600 at temperature where Sc(O2) ≈ 600", {
+  # Sc(O2) ≈ 600 when temp ≈ 17.55°C
+  k600 <- 2.0
+  temp_approx <- 17.55
+  Sc <- getSchmidt(temp_approx, "O2")
+  kGAS <- k600.2.kGAS.base(k600, temp_approx, "O2")
+  # kGAS should be close to k600 (within ~1%)
+  expect_equal(kGAS, k600, tolerance = 0.01)
+})
+
+test_that("k600.2.kGAS.base: larger gas molecules (higher Sc) yield smaller kGAS", {
+  # CO2 has higher Sc than O2 at same temperature, so kGAS should be smaller
+  k600 <- 2.0
+  kO2  <- k600.2.kGAS.base(k600, 20, "O2")
+  kCO2 <- k600.2.kGAS.base(k600, 20, "CO2")
+  expect_true(kCO2 < kO2)
+})
+
+test_that("k600.2.kGAS.base handles vector inputs", {
+  k600  <- c(1, 2, 3)
+  temps <- c(10, 15, 20)
+  result <- k600.2.kGAS.base(k600, temps, "O2")
+  expected <- k600 * (getSchmidt(temps, "O2") / 600)^(-0.5)
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("k600.2.kGAS wrapper returns data.frame with k.gas column", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00", "2020-01-01 01:00"), tz = "UTC"),
+    k600 = c(1.0, 2.0),
+    wtr  = c(15, 20))
+  result <- k600.2.kGAS(ts, gas = "O2")
+  expect_s3_class(result, "data.frame")
+  expect_true("k.gas" %in% names(result))
+  expect_true("datetime" %in% names(result))
+  expect_equal(result$k.gas,
+               k600.2.kGAS.base(ts$k600, ts$wtr, "O2"),
+               tolerance = 1e-10)
+})
diff --git a/tests/testthat/test_k_models.R b/tests/testthat/test_k_models.R
new file mode 100644
index 0000000..182376b
--- /dev/null
+++ b/tests/testthat/test_k_models.R
@@ -0,0 +1,117 @@
+## ---- k.cole ----
+
+test_that("k.cole.base known value at U10=4 m/s", {
+  # k600 = (2.07 + 0.215 * 4^1.7) * 24/100 in m/day
+  expected <- (2.07 + 0.215 * 4^1.7) * 24 / 100
+  result <- k.cole.base(4)
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("k.cole.base returns positive k600 at zero wind (intercept)", {
+  result <- k.cole.base(0)
+  expected <- 2.07 * 24 / 100
+  expect_equal(result, expected, tolerance = 1e-10)
+  expect_true(result > 0)
+})
+
+test_that("k.cole.base k600 increases with wind speed", {
+  expect_true(k.cole.base(5) > k.cole.base(3))
+  expect_true(k.cole.base(3) > k.cole.base(1))
+})
+
+test_that("k.cole wrapper returns data.frame with datetime and k600 columns", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00", "2020-01-01 01:00"), tz = "UTC"),
+    wnd = c(3, 5))
+  result <- k.cole(ts)
+  expect_s3_class(result, "data.frame")
+  expect_true(all(c("datetime", "k600") %in% names(result)))
+  expect_equal(result$k600, k.cole.base(ts$wnd))
+})
+
+test_that("k.cole errors when wnd column is missing", {
+  ts_bad <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wind = 5)
+  expect_error(k.cole(ts_bad))
+})
+
+## ---- k.crusius ----
+
+test_that("k.crusius.base power method returns numeric vector", {
+  result <- k.crusius.base(c(2, 4, 6), method = "power")
+  expect_type(result, "double")
+  expect_length(result, 3)
+  expect_true(all(is.finite(result)))
+})
+
+test_that("k.crusius.base constant method: all values below threshold equal same k600", {
+  # Below 3.7 m/s threshold → constant 1 cm/h = 0.24 m/day
+  result <- k.crusius.base(c(1, 2, 3), method = "constant")
+  expect_equal(result, rep(1 * 24 / 100, 3))
+})
+
+test_that("k.crusius.base constant method: value at threshold uses second formula", {
+  # At exactly 3.7, ifelse(3.7 < 3.7, ...) is FALSE → uses 5.14*3.7 - 17.9
+  expected <- (5.14 * 3.7 - 17.9) * 24 / 100
+  result <- k.crusius.base(3.7, method = "constant")
+  expect_equal(result, expected, tolerance = 1e-10)
+  expect_true(result > 1 * 24 / 100)  # higher than the below-threshold value
+})
+
+test_that("k.crusius.base bilinear method returns numeric vector", {
+  result <- k.crusius.base(c(2, 5), method = "bilinear")
+  expect_type(result, "double")
+  expect_length(result, 2)
+  # Below threshold
+  expect_equal(result[1], 0.72 * 2 * 24 / 100, tolerance = 1e-10)
+  # Above threshold
+  expect_equal(result[2], (4.33 * 5 - 13.3) * 24 / 100, tolerance = 1e-10)
+})
+
+test_that("k.crusius.base errors on invalid method", {
+  expect_error(k.crusius.base(5, method = "invalid"))
+})
+
+test_that("k.crusius wrapper returns data.frame with datetime and k600", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00", "2020-01-01 01:00"), tz = "UTC"),
+    wnd = c(3, 5))
+  result <- k.crusius(ts, method = "power")
+  expect_s3_class(result, "data.frame")
+  expect_true(all(c("datetime", "k600") %in% names(result)))
+  expect_equal(result$k600, k.crusius.base(ts$wnd, method = "power"))
+})
+
+## ---- k.vachon ----
+
+test_that("k.vachon.base known value at U10=3, lake.area=1e6 m2", {
+  # At 1e6 m^2, log10(1e6/1e6) = 0, so third term drops out
+  # k600 = (2.51 + 1.48*3 + 0) * 24/100
+  expected <- (2.51 + 1.48 * 3 + 0.39 * 3 * log10(1e6 / 1e6)) * 24 / 100
+  result <- k.vachon.base(3, 1e6)
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("k.vachon.base k600 increases with lake area", {
+  # log10(area/1e6) > 0 for area > 1e6 → adds to k600
+  expect_true(k.vachon.base(3, 1e7) > k.vachon.base(3, 1e6))
+  expect_true(k.vachon.base(3, 1e6) > k.vachon.base(3, 1e5))
+})
+
+test_that("k.vachon wrapper returns data.frame with datetime and k600", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00", "2020-01-01 01:00"), tz = "UTC"),
+    wnd = c(3, 5))
+  result <- k.vachon(ts, lake.area = 1e6)
+  expect_s3_class(result, "data.frame")
+  expect_true(all(c("datetime", "k600") %in% names(result)))
+  expect_equal(result$k600, k.vachon.base(ts$wnd, 1e6))
+})
+
+test_that("k.vachon errors when wnd column is missing", {
+  ts_bad <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wind = 5)
+  expect_error(k.vachon(ts_bad, lake.area = 1e6))
+})
diff --git a/tests/testthat/test_metab.R b/tests/testthat/test_metab.R
new file mode 100644
index 0000000..00d5146
--- /dev/null
+++ b/tests/testthat/test_metab.R
@@ -0,0 +1,102 @@
+# Shared synthetic data for metabolism tests:
+# 144 observations (10-min intervals) over one 24-hour period.
+make_metab_inputs <- function() {
+  n <- 144
+  start_time <- as.POSIXct("2020-07-15 00:00:00", tz = "UTC")
+  datetime   <- start_time + seq(0, (n - 1) * 600, by = 600)
+
+  wtr    <- rep(20, n)
+  do.sat <- o2.at.sat.base(wtr)
+  do.obs <- rep(8.5, n)   # slightly below saturation, no NAs
+  k.gas  <- rep(0.2, n)
+  z.mix  <- rep(2, n)
+
+  # Sinusoidal irr: 0 at night, peak ~1 at midday (hours 6–18 daytime)
+  hour_frac <- (seq_len(n) - 1) / 6   # 0 to 23.833 hours
+  irr <- pmax(0, sin(pi * (hour_frac - 6) / 12))
+
+  list(datetime = datetime, wtr = wtr, do.obs = do.obs, do.sat = do.sat,
+       k.gas = k.gas, z.mix = z.mix, irr = irr)
+}
+
+## ---- metab.ols ----
+
+test_that("metab.ols returns list with metab containing GPP, R, NEP", {
+  d <- make_metab_inputs()
+  result <- metab.ols(d$do.obs, d$do.sat, d$k.gas, d$z.mix, d$irr, d$wtr,
+                      datetime = d$datetime)
+  expect_type(result, "list")
+  expect_true("metab" %in% names(result))
+  metab <- result$metab
+  expect_true(all(c("GPP", "R", "NEP") %in% names(metab)))
+  expect_type(metab$GPP, "double")
+  expect_type(metab$R,   "double")
+  expect_type(metab$NEP, "double")
+  expect_true(is.finite(metab$GPP))
+  expect_true(is.finite(metab$R))
+  expect_true(is.finite(metab$NEP))
+})
+
+test_that("metab.ols NEP equals GPP + R", {
+  d <- make_metab_inputs()
+  result <- metab.ols(d$do.obs, d$do.sat, d$k.gas, d$z.mix, d$irr, d$wtr,
+                      datetime = d$datetime)
+  metab <- result$metab
+  expect_equal(metab$NEP, metab$GPP + metab$R, tolerance = 1e-10)
+})
+
+test_that("metab.ols errors when z.mix contains zero", {
+  d <- make_metab_inputs()
+  d$z.mix[1] <- 0
+  expect_error(metab.ols(d$do.obs, d$do.sat, d$k.gas, d$z.mix, d$irr, d$wtr))
+})
+
+test_that("metab.ols errors when inputs contain NA", {
+  d <- make_metab_inputs()
+  d$do.obs[10] <- NA
+  expect_error(metab.ols(d$do.obs, d$do.sat, d$k.gas, d$z.mix, d$irr, d$wtr))
+})
+
+## ---- metab.bookkeep ----
+
+test_that("metab.bookkeep returns data.frame with GPP, R, NEP columns", {
+  d <- make_metab_inputs()
+  # bookkeep requires integer 0/1 irr
+  irr_int <- as.integer(d$irr > 0)
+  result <- metab.bookkeep(d$do.obs, d$do.sat, d$k.gas, d$z.mix, irr_int,
+                            datetime = d$datetime)
+  expect_s3_class(result, "data.frame")
+  expect_true(all(c("GPP", "R", "NEP") %in% names(result)))
+  expect_type(result$GPP, "double")
+  expect_type(result$R,   "double")
+  expect_type(result$NEP, "double")
+})
+
+test_that("metab.bookkeep errors when irr is not integer 0/1 and datetime/lat absent", {
+  d <- make_metab_inputs()
+  # Continuous irr without datetime+lake.lat should error; suppress the
+  # "datetime not found" warning that is emitted before the stop().
+  expect_error(
+    suppressWarnings(
+      metab.bookkeep(d$do.obs, d$do.sat, d$k.gas, d$z.mix, d$irr)
+    )
+  )
+})
+
+test_that("metab.bookkeep R is non-positive for plausible night respiration", {
+  d <- make_metab_inputs()
+  irr_int <- as.integer(d$irr > 0)
+  result <- metab.bookkeep(d$do.obs, d$do.sat, d$k.gas, d$z.mix, irr_int,
+                            datetime = d$datetime)
+  # With constant do.obs slightly below sat, gas flux drives slight uptake.
+  # R is estimated from nighttime, should be <= 0 or at least numeric.
+  expect_true(is.numeric(result$R))
+})
+
+## ---- metab via file path (integration check) ----
+
+test_that("sparkling.doobs file exists and can be located", {
+  path <- system.file("extdata", "sparkling.doobs", package = "LakeMetabolizer")
+  expect_true(nchar(path) > 0)
+  expect_true(file.exists(path))
+})
diff --git a/tests/testthat/test_o2.at.sat.R b/tests/testthat/test_o2.at.sat.R
new file mode 100644
index 0000000..3a8defc
--- /dev/null
+++ b/tests/testthat/test_o2.at.sat.R
@@ -0,0 +1,55 @@
+test_that("o2.at.sat.base returns known value at 20C sea level freshwater (garcia-benson)", {
+  # Garcia-Benson at 20°C, sea level, freshwater ≈ 9.08 mg/L
+  result <- o2.at.sat.base(temp = 20, altitude = 0, salinity = 0,
+                            model = "garcia-benson")
+  expect_equal(result, 9.08, tolerance = 0.05)
+})
+
+test_that("o2.at.sat.base all four models return numeric values at 20C", {
+  models <- c("garcia-benson", "garcia", "weiss", "benson")
+  for (m in models) {
+    result <- o2.at.sat.base(temp = 20, altitude = 0, salinity = 0, model = m)
+    expect_type(result, "double")
+    expect_true(is.finite(result))
+  }
+})
+
+test_that("o2.at.sat.base saturation decreases with altitude", {
+  sat_low  <- o2.at.sat.base(temp = 20, altitude = 0)
+  sat_high <- o2.at.sat.base(temp = 20, altitude = 1000)
+  expect_true(sat_high < sat_low)
+})
+
+test_that("o2.at.sat.base saturation decreases with salinity", {
+  sat_fresh <- o2.at.sat.base(temp = 20, salinity = 0)
+  sat_saline <- o2.at.sat.base(temp = 20, salinity = 35)
+  expect_true(sat_saline < sat_fresh)
+})
+
+test_that("o2.at.sat.base saturation decreases with temperature", {
+  sat_cold <- o2.at.sat.base(temp = 5)
+  sat_warm <- o2.at.sat.base(temp = 25)
+  expect_true(sat_warm < sat_cold)
+})
+
+test_that("o2.at.sat.base stops on unrecognized model", {
+  expect_error(o2.at.sat.base(temp = 20, model = "unknown_model"))
+})
+
+test_that("o2.at.sat wrapper returns data.frame with do.sat column", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00:00", "2020-01-01 01:00:00"),
+                           tz = "UTC"),
+    wtr = c(10, 15))
+  result <- o2.at.sat(ts)
+  expect_s3_class(result, "data.frame")
+  expect_true("do.sat" %in% names(result))
+  expect_true("datetime" %in% names(result))
+  expect_equal(nrow(result), 2)
+  expect_equal(result$do.sat, o2.at.sat.base(ts$wtr))
+})
+
+test_that("benson model warns when salinity is non-zero", {
+  expect_warning(o2.at.sat.base(temp = 20, salinity = 10, model = "benson"),
+                 regexp = "Benson model does not currently include salinity")
+})
diff --git a/tests/testthat/test_par_conversions.R b/tests/testthat/test_par_conversions.R
new file mode 100644
index 0000000..3e068e0
--- /dev/null
+++ b/tests/testthat/test_par_conversions.R
@@ -0,0 +1,55 @@
+## ---- par.to.sw ----
+
+test_that("par.to.sw.base known value: 100 PAR -> 47.3 W/m2", {
+  expect_equal(par.to.sw.base(100), 47.3)
+})
+
+test_that("par.to.sw.base is linear in input", {
+  expect_equal(par.to.sw.base(200), 2 * par.to.sw.base(100))
+  expect_equal(par.to.sw.base(0),   0)
+})
+
+test_that("par.to.sw.base handles vector input", {
+  par_vals <- c(100, 200, 500)
+  expected <- par_vals * 0.473
+  expect_equal(par.to.sw.base(par_vals), expected)
+})
+
+test_that("par.to.sw wrapper returns data.frame with sw column", {
+  df <- data.frame(par = c(100, 500, 1000))
+  result <- par.to.sw(df)
+  expect_s3_class(result, "data.frame")
+  expect_true("sw" %in% names(result))
+  expect_false("par" %in% names(result))
+  expect_equal(result$sw, par.to.sw.base(c(100, 500, 1000)))
+})
+
+## ---- sw.to.par ----
+
+test_that("sw.to.par.base known value: 100 W/m2 -> 211.4 PAR", {
+  expect_equal(sw.to.par.base(100), 211.4)
+})
+
+test_that("sw.to.par.base is linear in input", {
+  expect_equal(sw.to.par.base(200), 2 * sw.to.par.base(100))
+  expect_equal(sw.to.par.base(0),   0)
+})
+
+test_that("sw.to.par wrapper returns data.frame with par column", {
+  df <- data.frame(sw = c(100, 500))
+  result <- sw.to.par(df)
+  expect_s3_class(result, "data.frame")
+  expect_true("par" %in% names(result))
+  expect_false("sw" %in% names(result))
+  expect_equal(result$par, sw.to.par.base(c(100, 500)))
+})
+
+## ---- roundtrip ----
+
+test_that("par -> sw -> par roundtrip is approximate (within 0.1%)", {
+  par_in  <- 100
+  sw_val  <- par.to.sw.base(par_in)
+  par_out <- sw.to.par.base(sw_val)
+  # 0.473 * 2.114 ≈ 0.99992, so roundtrip is not exact
+  expect_equal(par_out, par_in, tolerance = 1e-3)
+})
diff --git a/tests/testthat/test_schmidt.R b/tests/testthat/test_schmidt.R
new file mode 100644
index 0000000..66f8338
--- /dev/null
+++ b/tests/testthat/test_schmidt.R
@@ -0,0 +1,54 @@
+test_that("getSchmidt returns known value for O2 at 20C", {
+  # Sc = A + B*T + C*T^2 + D*T^3; O2 coefficients: 1568, -86.04, 2.142, -0.0216
+  expected <- 1568 + (-86.04) * 20 + 2.142 * 20^2 + (-0.0216) * 20^3
+  result <- getSchmidt(20, "O2")
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("getSchmidt returns numeric for all supported gases", {
+  gases <- c("He", "O2", "CO2", "CH4", "SF6", "N2O", "Ar", "N2")
+  for (g in gases) {
+    result <- getSchmidt(20, g)
+    expect_type(result, "double")
+    expect_true(is.finite(result))
+  }
+})
+
+test_that("getSchmidt errors on unrecognized gas", {
+  expect_error(getSchmidt(20, "Xe"))
+})
+
+test_that("getSchmidt warns when temperature is outside 4-35C range", {
+  expect_warning(getSchmidt(2, "O2"),  regexp = "temperature out of range")
+  expect_warning(getSchmidt(40, "O2"), regexp = "temperature out of range")
+})
+
+test_that("getSchmidt Schmidt number decreases as temperature increases (O2)", {
+  sc_low  <- getSchmidt(5,  "O2")
+  sc_mid  <- getSchmidt(15, "O2")
+  sc_high <- getSchmidt(30, "O2")
+  expect_true(sc_mid < sc_low)
+  expect_true(sc_high < sc_mid)
+})
+
+test_that("getSchmidt larger molecules have higher Schmidt numbers than smaller ones at same T", {
+  # SF6 (heavy) should have higher Sc than He (light) at same temperature
+  sc_He  <- getSchmidt(20, "He")
+  sc_SF6 <- getSchmidt(20, "SF6")
+  expect_true(sc_SF6 > sc_He)
+})
+
+test_that("getSchmidt handles vector temperature input", {
+  temps <- c(10, 20, 30)
+  result <- getSchmidt(temps, "O2")
+  expect_length(result, 3)
+  expect_type(result, "double")
+})
+
+test_that("getSchmidt handles NA in temperature vector without error", {
+  temps <- c(10, NA, 20)
+  # Should warn if non-NA values are out of range, but not error
+  result <- getSchmidt(temps, "O2")
+  expect_length(result, 3)
+  expect_true(is.na(result[2]))
+})
diff --git a/tests/testthat/test_sun.rise.set.R b/tests/testthat/test_sun.rise.set.R
index 13588a2..bb69ea1 100644
--- a/tests/testthat/test_sun.rise.set.R
+++ b/tests/testthat/test_sun.rise.set.R
@@ -1,22 +1,47 @@
-context("sun.rise.set()")
-library(LakeMetabolizer)
-
-testthat::test_that("Return object class and format", {
-  
-  # Run sun.rise.set() with test data
+test_that("sun.rise.set returns correct class and column names", {
   dates <- as.POSIXlt(
     c("2020-08-22", "2020-08-23", "2020-08-24"),
     tz = "America/Chicago")
   r <- sun.rise.set(dates, lat = 41.8781)
-  
-  # Returned object is a data frame
-  expect_equal(class(r), "data.frame")
-  
-  # Index order is preserved from original implmentation
-  expect_equal("sunrise", colnames(r)[1])
-  expect_equal("sunset", colnames(r)[2])
-  
-  # Input time zone attribute is returned
-  expect_true(attr(r$sunrise, "tzone") == "America/Chicago")
-  
+
+  expect_s3_class(r, "data.frame")
+  expect_equal(colnames(r)[1], "sunrise")
+  expect_equal(colnames(r)[2], "sunset")
+  expect_equal(attr(r$sunrise, "tzone"), "America/Chicago")
+})
+
+test_that("sun.rise.set sunrise is before sunset at mid-latitudes in summer", {
+  date <- as.POSIXct("2020-07-15 00:00:00", tz = "UTC")
+  r <- sun.rise.set(date, lat = 45)
+  expect_true(r$sunrise < r$sunset)
+})
+
+test_that("sun.rise.set returns NA at high latitude polar day", {
+  # At lat=89 in summer, sun does not set: omegaInput < -1 → NA
+  summer_date <- as.POSIXct("2020-06-21", tz = "UTC")
+  r <- sun.rise.set(summer_date, lat = 89)
+  expect_true(is.na(r$sunrise))
+  expect_true(is.na(r$sunset))
+})
+
+test_that("is.day and is.night are complementary boolean vectors", {
+  datetimes <- as.POSIXct(
+    c("2020-07-15 06:00:00", "2020-07-15 14:00:00", "2020-07-15 22:00:00"),
+    tz = "UTC")
+  lat <- 45
+  day_vals   <- is.day(datetimes, lat)
+  night_vals <- is.night(datetimes, lat)
+
+  expect_type(day_vals, "logical")
+  expect_type(night_vals, "logical")
+  expect_length(day_vals, 3)
+  expect_equal(day_vals, !night_vals)
+})
+
+test_that("is.day returns TRUE at midday and FALSE at midnight", {
+  noon     <- as.POSIXct("2020-07-15 12:00:00", tz = "UTC")
+  midnight <- as.POSIXct("2020-07-15 00:00:00", tz = "UTC")
+  lat <- 45
+  expect_true(is.day(noon, lat))
+  expect_false(is.day(midnight, lat))
 })
diff --git a/tests/testthat/test_wind.scale.R b/tests/testthat/test_wind.scale.R
new file mode 100644
index 0000000..7966287
--- /dev/null
+++ b/tests/testthat/test_wind.scale.R
@@ -0,0 +1,47 @@
+test_that("wind.scale.base returns correct value: wnd=5 at height=2", {
+  # U10 = wnd * (10/wnd.z)^0.15 = 5 * 5^0.15
+  expected <- 5 * (10 / 2)^0.15
+  result <- wind.scale.base(5, 2)
+  expect_equal(result, expected, tolerance = 1e-10)
+})
+
+test_that("wind.scale.base identity: wind at 10m is unchanged", {
+  result <- wind.scale.base(5, 10)
+  expect_equal(result, 5)
+})
+
+test_that("wind.scale.base scales up wind from height below 10m", {
+  result <- wind.scale.base(5, 2)
+  expect_true(result > 5)
+})
+
+test_that("wind.scale.base scales down wind from height above 10m", {
+  result <- wind.scale.base(5, 20)
+  expect_true(result < 5)
+})
+
+test_that("wind.scale.base handles vector input", {
+  wnd <- c(3, 5, 7)
+  result <- wind.scale.base(wnd, 2)
+  expected <- wnd * (10 / 2)^0.15
+  expect_equal(result, expected)
+})
+
+test_that("wind.scale wrapper returns data.frame with wnd_10 column", {
+  ts <- data.frame(
+    datetime = as.POSIXct(c("2020-01-01 00:00:00", "2020-01-01 01:00:00"),
+                           tz = "UTC"),
+    wnd = c(3, 5))
+  result <- wind.scale(ts, wnd.z = 2)
+  expect_s3_class(result, "data.frame")
+  expect_true("wnd_10" %in% names(result))
+  expect_true("datetime" %in% names(result))
+  expect_equal(result$wnd_10, wind.scale.base(ts$wnd, 2))
+})
+
+test_that("wind.scale errors when wnd column is missing", {
+  ts_bad <- data.frame(
+    datetime = as.POSIXct("2020-01-01", tz = "UTC"),
+    wind_speed = 5)
+  expect_error(wind.scale(ts_bad, wnd.z = 2))
+})