diff --git a/DESCRIPTION b/DESCRIPTION
index fc9bdf9..de76b4e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -31,7 +31,8 @@ Depends:
     R (>= 4.1)
 Imports:
     terra,
-    sf
+    sf,
+    tibble
 Suggests:
     bookdown,
     whitebox,
diff --git a/NAMESPACE b/NAMESPACE
index a1696aa..7cd9c3f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,8 +8,10 @@ export(fl_flood_surface)
 export(fl_flood_trim)
 export(fl_mask)
 export(fl_mask_distance)
+export(fl_params)
 export(fl_patch_conn)
 export(fl_patch_rm)
+export(fl_scenarios)
 export(fl_stream_rasterize)
 export(fl_valley_confine)
 export(fl_valley_poly)
diff --git a/R/fl_params.R b/R/fl_params.R
new file mode 100644
index 0000000..366c9e9
--- /dev/null
+++ b/R/fl_params.R
@@ -0,0 +1,29 @@
+#' Load VCA parameter legend
+#'
+#' Returns a tibble of Valley Confinement Algorithm parameters with units,
+#' defaults, literature sources, and descriptions. The bundled CSV documents
+#' every tuning parameter in [fl_valley_confine()] with DEM resolution
+#' guidance from Nagel et al. (2014) and Hall et al. (2007).
+#'
+#' @param path Character. Path to a custom parameter CSV. When `NULL`
+#'   (default), loads the bundled `inst/extdata/flood_params.csv`.
+#'
+#' @return A tibble with columns: `parameter`, `unit`, `default`, `source`,
+#'   `citation_keys`, `effect`, `description`.
+#'
+#' @examples
+#' params <- fl_params()
+#' params[, c("parameter", "unit", "default")]
+#'
+#' @export
+fl_params <- function(path = NULL) {
+  if (is.null(path)) {
+    path <- system.file("extdata", "flood_params.csv", package = "flooded",
+                        mustWork = TRUE)
+  }
+  if (!file.exists(path)) {
+    stop("File not found: ", path, call. = FALSE)
+  }
+  utils::read.csv(path, stringsAsFactors = FALSE) |>
+    tibble::as_tibble()
+}
diff --git a/R/fl_scenarios.R b/R/fl_scenarios.R
new file mode 100644
index 0000000..350cdd2
--- /dev/null
+++ b/R/fl_scenarios.R
@@ -0,0 +1,52 @@
+#' Load flood factor scenarios
+#'
+#' Returns a tibble of pre-defined flood factor scenarios for
+#' [fl_valley_confine()]. Each row is a complete parameter set that can be
+#' passed to the function. The bundled CSV includes three scenarios spanning
+#' active channel margin (ff=2) to full valley bottom (ff=6).
+#'
+#' The `flood_factor` is a DEM compensation parameter, not an ecological
+#' threshold — no paper maps specific ff values to ecological processes.
+#' The scenario descriptions are an interpretive overlay based on where
+#' different ff values fall relative to field-validated studies:
+#'
+#' - **ff=2**: Rosgen flood-prone width, ~50-yr flood stage
+#' - **ff=3-4**: Historical floodplain (Hall et al. 2007 validated ff=3 on
+#'   10m DEM against 213 field sites; ff=4 compensates for 25m DEM smoothing)
+#' - **ff=5-7**: Valley bottom including terraces (Nagel et al. 2014)
+#'
+#' DEM resolution matters: coarser DEMs need larger ff to compensate for
+#' smoothed valley floors. At 1m lidar, ff=2-3 may suffice.
+#'
+#' The `run` column allows consuming projects to document all scenarios but
+#' only execute selected ones (e.g., `dplyr::filter(scenarios, run)`).
+#'
+#' @param path Character. Path to a custom scenarios CSV. When `NULL`
+#'   (default), loads the bundled `inst/extdata/flood_scenarios.csv`.
+#'   Users can copy the default CSV, add rows for site-level work
+#'   (e.g., smaller thresholds for 1m lidar), and pass the custom path.
+#'
+#' @return A tibble with columns: `scenario_id`, `flood_factor`,
+#'   `slope_threshold`, `max_width`, `cost_threshold`, `size_threshold`,
+#'   `hole_threshold`, `run`, `description`, `ecological_process`,
+#'   `citation_keys`.
+#'
+#' @examples
+#' scenarios <- fl_scenarios()
+#' scenarios[, c("scenario_id", "flood_factor", "description")]
+#'
+#' # Filter to scenarios marked for execution
+#' to_run <- scenarios[scenarios$run, ]
+#'
+#' @export
+fl_scenarios <- function(path = NULL) {
+  if (is.null(path)) {
+    path <- system.file("extdata", "flood_scenarios.csv", package = "flooded",
+                        mustWork = TRUE)
+  }
+  if (!file.exists(path)) {
+    stop("File not found: ", path, call. = FALSE)
+  }
+  utils::read.csv(path, stringsAsFactors = FALSE) |>
+    tibble::as_tibble()
+}
diff --git a/inst/extdata/flood_params.csv b/inst/extdata/flood_params.csv
index 947092b..afffb37 100644
--- a/inst/extdata/flood_params.csv
+++ b/inst/extdata/flood_params.csv
@@ -1,7 +1,7 @@
-parameter,unit,default,source,citation_keys,description
-flood_factor,dimensionless,6,"Nagel et al. 2014; Hall et al. 2007",nagel_etal2014LandscapeScale;hall_etal2007Predictingriver,"Bankfull depth multiplier — scales predicted bankfull depth to set flood height above channel. Rosgen ff=2 approximates 50-yr flood. @hall_etal2007Predictingriver found ff=3 best fit for historical floodplain on 10m DEM (213 field sites). @nagel_etal2014LandscapeScale default ff=5 for 10m DEM and ff=7 for 30m DEM — these map valley bottom (wider than floodplain). Higher values compensate for coarser DEM vertical resolution. At 1m LiDAR ff=2-3 may suffice (minimal smoothing). Critical distinction: ff=3-4 approximates functional floodplain; ff=5-7 maps full valley bottom including terraces."
-slope_threshold,percent,9,Nagel et al. 2014,nagel_etal2014LandscapeScale,"Maximum percent slope for valley floor cells. @nagel_etal2014LandscapeScale: derived empirically from field surveys in central Idaho — slopes <9% in DEM likely correspond to unconfined valleys. Note: 9% slope (~5.1 degrees) appears high but DEM-derived slopes are systematically steeper than true ground values due to rounding of valley bottoms. Calibrated on 30m DEMs. At 1m LiDAR where slopes are more accurate a lower threshold may be appropriate."
-max_width,metres,2000,Nagel et al. 2014 (modified),nagel_etal2014LandscapeScale,"Maximum valley width orthogonal to channel. @nagel_etal2014LandscapeScale called this an arbitrary measure and used 500-1000m. Safety cap preventing runaway flood spread in flat terrain — in most valleys flood_factor and slope_threshold constrain extent before this kicks in. Increase for wide mainstem valleys. Not sensitive to DEM resolution."
-cost_threshold,dimensionless,2500,Nagel et al. 2014,nagel_etal2014LandscapeScale,"Accumulated cost distance limit (cell_size x cell_slope summed outward from stream). @nagel_etal2014LandscapeScale: has no physical meaning — empirical rule setting initial processing domain. Intentionally overestimates valley extent for subsequent refinement by flood_factor and slope_threshold. Scales with DEM resolution: at 1m the same physical distance accumulates ~25x less cost than at 25m. May need adjustment for high-resolution DEMs."
-size_threshold,m2,5000,Nagel et al. 2014 (modified),nagel_etal2014LandscapeScale,"Minimum valley patch area — removes noise patches smaller than this. @nagel_etal2014LandscapeScale used 10000 m2 (1 ha). flooded default 5000 m2 (0.5 ha). At 25m DEM = 8 pixels minimum. At 1m LiDAR = 5000 pixels — preserves very small features. For 1m work reduce to 100-500 m2 to filter noise while keeping side channels and small floodplain features."
-hole_threshold,m2,2500,flooded-specific,NA,"Maximum hole area to fill in floodplain polygon. Bridges gaps from roads or infrastructure. flooded-specific — not from @nagel_etal2014LandscapeScale. At 25m DEM = 4 pixels (fills small DEM artifact gaps). At 1m LiDAR = 2500 pixels (~50m x 50m) — large enough to erase roads and berms. For site-level restoration design where infrastructure visibility matters reduce to 25-100 m2."
+parameter,unit,default,source,citation_keys,effect,description
+flood_factor,dimensionless,6,"Nagel et al. 2014; Hall et al. 2007",nagel_etal2014LandscapeScale;hall_etal2007Predictingriver,Higher = deeper flood; more floodplain,"Bankfull depth multiplier — scales predicted bankfull depth to set flood height above channel. Rosgen ff=2 approximates 50-yr flood. @hall_etal2007Predictingriver found ff=3 best fit for historical floodplain on 10m DEM (213 field sites). @nagel_etal2014LandscapeScale default ff=5 for 10m DEM and ff=7 for 30m DEM — these map valley bottom (wider than floodplain). Higher values compensate for coarser DEM vertical resolution. At 1m LiDAR ff=2-3 may suffice (minimal smoothing). Critical distinction: ff=3-4 approximates functional floodplain; ff=5-7 maps full valley bottom including terraces."
+slope_threshold,percent,9,Nagel et al. 2014,nagel_etal2014LandscapeScale,Higher = more valley floor included,"Maximum percent slope for valley floor cells. @nagel_etal2014LandscapeScale: derived empirically from field surveys in central Idaho — slopes <9% in DEM likely correspond to unconfined valleys. Note: 9% slope (~5.1 degrees) appears high but DEM-derived slopes are systematically steeper than true ground values due to rounding of valley bottoms. Calibrated on 30m DEMs. At 1m LiDAR where slopes are more accurate a lower threshold may be appropriate."
+max_width,metres,2000,Nagel et al. 2014 (modified),nagel_etal2014LandscapeScale,Analysis corridor width,"Maximum valley width orthogonal to channel. @nagel_etal2014LandscapeScale called this an arbitrary measure and used 500-1000m. Safety cap preventing runaway flood spread in flat terrain — in most valleys flood_factor and slope_threshold constrain extent before this kicks in. Increase for wide mainstem valleys. Not sensitive to DEM resolution."
+cost_threshold,dimensionless,2500,Nagel et al. 2014,nagel_etal2014LandscapeScale,Higher = valley extends further up hillslopes,"Accumulated cost distance limit (cell_size x cell_slope summed outward from stream). @nagel_etal2014LandscaleScale: has no physical meaning — empirical rule setting initial processing domain. Intentionally overestimates valley extent for subsequent refinement by flood_factor and slope_threshold. Scales with DEM resolution: at 1m the same physical distance accumulates ~25x less cost than at 25m. May need adjustment for high-resolution DEMs."
+size_threshold,m2,5000,Nagel et al. 2014 (modified),nagel_etal2014LandscapeScale,Minimum valley patch area,"Minimum valley patch area — removes noise patches smaller than this. @nagel_etal2014LandscapeScale used 10000 m2 (1 ha). flooded default 5000 m2 (0.5 ha). At 25m DEM = 8 pixels minimum. At 1m LiDAR = 5000 pixels — preserves very small features. For 1m work reduce to 100-500 m2 to filter noise while keeping side channels and small floodplain features."
+hole_threshold,m2,2500,flooded-specific,NA,Maximum hole size to fill,"Maximum hole area to fill in floodplain polygon. Bridges gaps from roads or infrastructure. flooded-specific — not from @nagel_etal2014LandscapeScale. At 25m DEM = 4 pixels (fills small DEM artifact gaps). At 1m LiDAR = 2500 pixels (~50m x 50m) — large enough to erase roads and berms. For site-level restoration design where infrastructure visibility matters reduce to 25-100 m2."
diff --git a/inst/extdata/flood_scenarios.csv b/inst/extdata/flood_scenarios.csv
new file mode 100644
index 0000000..722b321
--- /dev/null
+++ b/inst/extdata/flood_scenarios.csv
@@ -0,0 +1,4 @@
+scenario_id,flood_factor,slope_threshold,max_width,cost_threshold,size_threshold,hole_threshold,run,description,ecological_process,citation_keys
+ff02,2,9,2000,2500,5000,2500,TRUE,Flood-prone width / active channel margin,Rosgen flood-prone width approximating 50-yr flood stage. Captures the zone of frequent inundation and active channel migration.,nagel_etal2014LandscapeScale
+ff04,4,9,2000,2500,5000,2500,TRUE,Functional floodplain,Historical floodplain extent where nutrient exchange and LWD recruitment occur. Hall et al. found ff=3 best fit on 10m DEM; ff=4 compensates for 25m TRIM vertical smoothing.,hall_etal2007Predictingriver;nagel_etal2014LandscapeScale
+ff06,6,9,2000,2500,5000,2500,TRUE,Valley bottom extent,Full depositional zone including terraces. Nagel et al. recommended ff=5-7 for valley bottom mapping. Wider than functional floodplain — includes areas not regularly influenced by high flows.,nagel_etal2014LandscapeScale
diff --git a/man/fl_params.Rd b/man/fl_params.Rd
new file mode 100644
index 0000000..bfc74b6
--- /dev/null
+++ b/man/fl_params.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fl_params.R
+\name{fl_params}
+\alias{fl_params}
+\title{Load VCA parameter legend}
+\usage{
+fl_params(path = NULL)
+}
+\arguments{
+\item{path}{Character. Path to a custom parameter CSV. When \code{NULL}
+(default), loads the bundled \code{inst/extdata/flood_params.csv}.}
+}
+\value{
+A tibble with columns: \code{parameter}, \code{unit}, \code{default}, \code{source},
+\code{citation_keys}, \code{effect}, \code{description}.
+}
+\description{
+Returns a tibble of Valley Confinement Algorithm parameters with units,
+defaults, literature sources, and descriptions. The bundled CSV documents
+every tuning parameter in \code{\link[=fl_valley_confine]{fl_valley_confine()}} with DEM resolution
+guidance from Nagel et al. (2014) and Hall et al. (2007).
+}
+\examples{
+params <- fl_params()
+params[, c("parameter", "unit", "default")]
+
+}
diff --git a/man/fl_scenarios.Rd b/man/fl_scenarios.Rd
new file mode 100644
index 0000000..a7c3b57
--- /dev/null
+++ b/man/fl_scenarios.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fl_scenarios.R
+\name{fl_scenarios}
+\alias{fl_scenarios}
+\title{Load flood factor scenarios}
+\usage{
+fl_scenarios(path = NULL)
+}
+\arguments{
+\item{path}{Character. Path to a custom scenarios CSV. When \code{NULL}
+(default), loads the bundled \code{inst/extdata/flood_scenarios.csv}.
+Users can copy the default CSV, add rows for site-level work
+(e.g., smaller thresholds for 1m lidar), and pass the custom path.}
+}
+\value{
+A tibble with columns: \code{scenario_id}, \code{flood_factor},
+\code{slope_threshold}, \code{max_width}, \code{cost_threshold}, \code{size_threshold},
+\code{hole_threshold}, \code{run}, \code{description}, \code{ecological_process},
+\code{citation_keys}.
+}
+\description{
+Returns a tibble of pre-defined flood factor scenarios for
+\code{\link[=fl_valley_confine]{fl_valley_confine()}}. Each row is a complete parameter set that can be
+passed to the function. The bundled CSV includes three scenarios spanning
+active channel margin (ff=2) to full valley bottom (ff=6).
+}
+\details{
+The \code{flood_factor} is a DEM compensation parameter, not an ecological
+threshold — no paper maps specific ff values to ecological processes.
+The scenario descriptions are an interpretive overlay based on where
+different ff values fall relative to field-validated studies:
+\itemize{
+\item \strong{ff=2}: Rosgen flood-prone width, ~50-yr flood stage
+\item \strong{ff=3-4}: Historical floodplain (Hall et al. 2007 validated ff=3 on
+10m DEM against 213 field sites; ff=4 compensates for 25m DEM smoothing)
+\item \strong{ff=5-7}: Valley bottom including terraces (Nagel et al. 2014)
+}
+
+DEM resolution matters: coarser DEMs need larger ff to compensate for
+smoothed valley floors. At 1m lidar, ff=2-3 may suffice.
+
+The \code{run} column allows consuming projects to document all scenarios but
+only execute selected ones (e.g., \code{dplyr::filter(scenarios, run)}).
+}
+\examples{
+scenarios <- fl_scenarios()
+scenarios[, c("scenario_id", "flood_factor", "description")]
+
+# Filter to scenarios marked for execution
+to_run <- scenarios[scenarios$run, ]
+
+}
diff --git a/planning/active/findings.md b/planning/active/findings.md
new file mode 100644
index 0000000..8b4385a
--- /dev/null
+++ b/planning/active/findings.md
@@ -0,0 +1,27 @@
+# Findings: fl_scenarios() and fl_params()
+
+## Research (pre-existing)
+
+All research already committed:
+- `inst/research/vca_parameter_rationale.md` — verified quotes from Nagel 2014, Hall 2007
+- `inst/extdata/flood_params.csv` — parameter legend with units, defaults, citations
+
+## Key facts for scenario design
+
+### flood_factor zones (from vca_parameter_rationale.md)
+- ff=2: Rosgen flood-prone width, ~50-yr flood (active channel margin)
+- ff=3: Hall 2007 best fit for historical floodplain on 10m DEM (213 field sites)
+- ff=4: Reasonable functional floodplain estimate on 25m DEM
+- ff=5-7: Nagel 2014 valley bottom (includes terraces, depositional areas)
+- ff=7: Nagel 2014 recommendation for 30m DEM
+
+### Defaults from flood_params.csv
+- slope_threshold: 9 (percent)
+- max_width: 2000 (metres)
+- cost_threshold: 2500 (dimensionless)
+- size_threshold: 5000 (m²)
+- hole_threshold: 2500 (m²)
+
+### Two independent width models (from flooded#25)
+- bcfishpass (Thorley et al. 2021): `exp(0.307) * (area * precip / 100000) ^ 0.458`
+- VCA/Hall 2007 (flooded internal): `(area ^ 0.280) * 0.196 * (precip ^ 0.355)`
diff --git a/planning/active/progress.md b/planning/active/progress.md
new file mode 100644
index 0000000..38fe28b
--- /dev/null
+++ b/planning/active/progress.md
@@ -0,0 +1,11 @@
+# Progress: fl_scenarios() and fl_params()
+
+## Session log
+
+### 2026-03-19
+- Read flooded#28 comments — updated scope and implementation prompt
+- Read inst/research/vca_parameter_rationale.md and inst/extdata/flood_params.csv
+- Read R/fl_valley_confine.R for current param defaults
+- Created branch `flood-scenarios` from main
+- Set up planning/active/ structure
+- Starting Phase 1: flood_scenarios.csv
diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md
new file mode 100644
index 0000000..78f3777
--- /dev/null
+++ b/planning/active/task_plan.md
@@ -0,0 +1,48 @@
+# Task: fl_scenarios() and fl_params() — CSV loaders for VCA parameters
+
+**Issue:** flooded#28
+**Branch:** flood-scenarios
+
+## Goal
+Add two exported functions (`fl_scenarios()`, `fl_params()`) and one new CSV
+(`flood_scenarios.csv`) to give users access to VCA parameter metadata and
+pre-defined flood factor scenarios. Close #28.
+
+## Phases
+
+### Phase 1: Create flood_scenarios.csv
+- [x] Create `inst/extdata/flood_scenarios.csv`
+- [x] Columns: scenario_id, flood_factor, slope_threshold, max_width, cost_threshold, size_threshold, hole_threshold, run, description, ecological_process, citation_keys
+- [x] Three rows: ff02 (active channel), ff04 (functional floodplain), ff06 (valley bottom)
+- [x] Non-flood_factor params use defaults from flood_params.csv
+
+### Phase 2: fl_params()
+- [x] Create `R/fl_params.R`
+- [x] Load flood_params.csv by default, accept custom path
+- [x] Return tibble
+- [x] roxygen with @export
+
+### Phase 3: fl_scenarios()
+- [x] Create `R/fl_scenarios.R`
+- [x] Load flood_scenarios.csv by default, accept custom path
+- [x] Return tibble
+- [x] roxygen with @export
+
+### Phase 4: Tests
+- [x] Create `tests/testthat/test-fl_params.R` (4 tests)
+- [x] Create `tests/testthat/test-fl_scenarios.R` (7 tests)
+- [x] Default load returns tibble with expected columns
+- [x] Custom path override works
+- [x] Invalid path errors cleanly
+- [x] Scenarios match flood_params.csv defaults for non-flood_factor params
+
+### Phase 5: Document and verify
+- [x] devtools::document() — NAMESPACE, fl_params.Rd, fl_scenarios.Rd
+- [x] devtools::test() — 154 pass (22 new)
+- [x] lintr::lint_package() — no new warnings
+- [x] Commit with Fixes #28 (16bf126)
+
+## Constraints
+- Do NOT modify fl_valley_confine()
+- Do NOT add vignette (separate issue)
+- Do NOT bump version
diff --git a/tests/testthat/test-fl_params.R b/tests/testthat/test-fl_params.R
new file mode 100644
index 0000000..25a2512
--- /dev/null
+++ b/tests/testthat/test-fl_params.R
@@ -0,0 +1,33 @@
+test_that("fl_params returns tibble with expected columns", {
+  params <- fl_params()
+  expect_s3_class(params, "tbl_df")
+  expected_cols <- c("parameter", "unit", "default", "source",
+                     "citation_keys", "effect", "description")
+  expect_true(all(expected_cols %in% names(params)))
+})
+
+test_that("fl_params has all fl_valley_confine parameters", {
+  params <- fl_params()
+  expected_params <- c("flood_factor", "slope_threshold", "max_width",
+                       "cost_threshold", "size_threshold", "hole_threshold")
+  expect_true(all(expected_params %in% params$parameter))
+})
+
+test_that("fl_params accepts custom path", {
+  # Write a minimal CSV
+  tmp <- tempfile(fileext = ".csv")
+  on.exit(unlink(tmp))
+  writeLines(
+    c("parameter,unit,default,source,citation_keys,description",
+      "test_param,m,10,test,NA,A test parameter"),
+    tmp
+  )
+  params <- fl_params(path = tmp)
+  expect_s3_class(params, "tbl_df")
+  expect_equal(nrow(params), 1L)
+  expect_equal(params$parameter, "test_param")
+})
+
+test_that("fl_params errors on invalid path", {
+  expect_error(fl_params(path = "/nonexistent/file.csv"), "not found")
+})
diff --git a/tests/testthat/test-fl_scenarios.R b/tests/testthat/test-fl_scenarios.R
new file mode 100644
index 0000000..7deff19
--- /dev/null
+++ b/tests/testthat/test-fl_scenarios.R
@@ -0,0 +1,59 @@
+test_that("fl_scenarios returns tibble with expected columns", {
+  scenarios <- fl_scenarios()
+  expect_s3_class(scenarios, "tbl_df")
+  expected_cols <- c("scenario_id", "flood_factor", "slope_threshold",
+                     "max_width", "cost_threshold", "size_threshold",
+                     "hole_threshold", "run", "description",
+                     "ecological_process", "citation_keys")
+  expect_true(all(expected_cols %in% names(scenarios)))
+})
+
+test_that("fl_scenarios has three default scenarios", {
+  scenarios <- fl_scenarios()
+  expect_equal(nrow(scenarios), 3L)
+  expect_equal(scenarios$scenario_id, c("ff02", "ff04", "ff06"))
+})
+
+test_that("fl_scenarios flood_factor values match scenario IDs", {
+  scenarios <- fl_scenarios()
+  expect_equal(scenarios$flood_factor, c(2, 4, 6))
+})
+
+test_that("fl_scenarios non-flood_factor params match defaults", {
+  scenarios <- fl_scenarios()
+  params <- fl_params()
+  for (p in c("slope_threshold", "max_width", "cost_threshold",
+              "size_threshold", "hole_threshold")) {
+    default_val <- as.numeric(params$default[params$parameter == p])
+    expect_true(all(scenarios[[p]] == default_val),
+                info = paste(p, "should match flood_params.csv default"))
+  }
+})
+
+test_that("fl_scenarios run column is logical", {
+  scenarios <- fl_scenarios()
+  expect_type(scenarios$run, "logical")
+})
+
+test_that("fl_scenarios accepts custom path", {
+  tmp <- tempfile(fileext = ".csv")
+  on.exit(unlink(tmp))
+  header <- paste0(
+    "scenario_id,flood_factor,slope_threshold,max_width,",
+    "cost_threshold,size_threshold,hole_threshold,run,",
+    "description,ecological_process,citation_keys"
+  )
+  writeLines(
+    c(header,
+      "custom,3,9,2000,2500,5000,2500,TRUE,Custom scenario,Test,NA"),
+    tmp
+  )
+  scenarios <- fl_scenarios(path = tmp)
+  expect_s3_class(scenarios, "tbl_df")
+  expect_equal(nrow(scenarios), 1L)
+  expect_equal(scenarios$scenario_id, "custom")
+})
+
+test_that("fl_scenarios errors on invalid path", {
+  expect_error(fl_scenarios(path = "/nonexistent/file.csv"), "not found")
+})
diff --git a/vignettes/valley-confinement.Rmd b/vignettes/valley-confinement.Rmd
index cbea836..476e1a2 100644
--- a/vignettes/valley-confinement.Rmd
+++ b/vignettes/valley-confinement.Rmd
@@ -184,11 +184,12 @@ and the resulting valley is significantly narrower.
 valleys <- fl_valley_confine(
   dem, streams,
   field = "upstream_area_ha",
-  slope_threshold = 9,
-  max_width = 2000,
-  cost_threshold = 2500,
-  flood_factor = 6,
-  precip = precip_r
+  slope = slope,            # pre-computed; derived from DEM if NULL
+  slope_threshold = 9,      # percent slope — cells steeper are excluded
+  max_width = 2000,         # metres — safety cap on valley width
+  cost_threshold = 2500,    # accumulated cost distance limit
+  flood_factor = 6,         # bankfull depth multiplier (valley bottom)
+  precip = precip_r         # mean annual precipitation (mm)
 )
 
 n_valley <- sum(values(valleys) == 1, na.rm = TRUE)
@@ -336,6 +337,77 @@ cat("With precip:   ", n_valley, "cells (",
     round(100 * n_valley / ncell(dem), 1), "%)\n")
 ```
 
+## Flood factor scenarios
+
+The `flood_factor` is a DEM compensation parameter — it scales predicted
+bankfull depth to set the flood height above the channel. Higher values
+compensate for coarser DEM vertical resolution. The ecological interpretation
+is an overlay: Hall et al. (2007) validated ff=3 against 213 field-measured
+floodplain widths on 10 m DEM; Nagel et al. (2014) recommended ff=5–7 for
+mapping valley bottoms (wider than functional floodplain).
+
+`fl_scenarios()` provides three pre-defined scenarios. `fl_params()` documents
+every tuning parameter with units, defaults, and literature sources.
+
+```{r scenarios}
+scenarios <- fl_scenarios()
+scenarios[, c("scenario_id", "flood_factor", "description")]
+```
+
+```{r params}
+params <- fl_params()
+params[, c("parameter", "unit", "default", "source")]
+```
+
+Run all three scenarios on the test data to see how flood factor controls
+the mapped extent:
+
+```{r run-scenarios}
+results <- list()
+for (i in seq_len(nrow(scenarios))) {
+  s <- scenarios[i, ]
+  results[[s$scenario_id]] <- fl_valley_confine(
+    dem, streams,
+    field = "upstream_area_ha",
+    slope = slope,
+    precip = precip_r,
+    flood_factor = s$flood_factor,
+    slope_threshold = s$slope_threshold,
+    max_width = s$max_width,
+    cost_threshold = s$cost_threshold,
+    size_threshold = s$size_threshold,
+    hole_threshold = s$hole_threshold,
+    waterbodies = waterbodies
+  )
+}
+
+cell_area <- prod(res(dem))
+for (id in names(results)) {
+  n <- sum(values(results[[id]]) == 1, na.rm = TRUE)
+  cat(sprintf("%-6s (ff=%s): %5.1f ha\n", id,
+              scenarios$flood_factor[scenarios$scenario_id == id],
+              n * cell_area / 1e4))
+}
+```
+
+```{r plot-scenarios, fig.width = 7, fig.height = 10, fig.cap = "Valley delineation at three flood factor values. ff=2 (active channel margin) captures the zone of frequent inundation. ff=4 (functional floodplain) maps the historical flood extent. ff=6 (valley bottom) includes terraces and depositional surfaces beyond the active floodplain. Blue lines = streams, orange outlines = waterbody polygons."}
+par(mfrow = c(3, 1), mar = c(2, 4, 2, 1))
+titles <- c(
+  ff02 = "ff=2: Active channel margin",
+  ff04 = "ff=4: Functional floodplain",
+  ff06 = "ff=6: Valley bottom"
+)
+for (id in names(results)) {
+  plot(results[[id]], col = c("grey90", "darkgreen"),
+       main = titles[id], legend = FALSE)
+  plot(st_geometry(streams), add = TRUE, col = "blue", lwd = 1)
+  if (nrow(waterbodies) > 0) {
+    plot(st_geometry(waterbodies), add = TRUE,
+         border = "darkorange", lwd = 1.2)
+  }
+}
+```
+
 ## Performance
 
 Several `terra` operations inside `fl_valley_confine()` support multi-threading
@@ -352,16 +424,20 @@ which are single-threaded.
 
 ## Summary
 
-Key tuning parameters:
+Key tuning parameters from `fl_params()`:
+
+```{r summary-table}
+params <- fl_params()
+knitr::kable(
+  params[, c("parameter", "default", "unit", "effect")],
+  col.names = c("Parameter", "Default", "Unit", "Effect")
+)
+```
+
+Additional arguments to `fl_valley_confine()` not in the parameter legend:
 
 | Parameter | Default | Effect |
 |-----------|---------|--------|
-| `slope_threshold` | 9% | Higher = more valley floor included |
-| `max_width` | 2000 m | Analysis corridor width |
-| `cost_threshold` | 2500 | Higher = valley extends further up hillslopes |
-| `flood_factor` | 6 | Higher = deeper flood, more floodplain |
 | `precip` | 1 | MAP in mm — critical for realistic flood depth |
 | `waterbodies` | NULL | sf polygons of lakes/wetlands to fill donut holes |
 | `channel_buffer` | auto | Buffer streams by channel_width (DEM correction) |
-| `size_threshold` | 5000 m² | Minimum valley patch area |
-| `hole_threshold` | 2500 m² | Maximum hole size to fill |