diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9b14f4ebe4e..83f2c3f61a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 ## Unreleased
 
 ### Added
+- Added `sensobol` based global sensitivity design and postprocessing support in `PEcAn.uncertainty`, including saved Sobol design and index outputs.
 - Added PEcAn.PEPRMT model, including a demo run with example data
 - Add `format_try_for_ma()` and `try_trait_mapping()` to `PEcAn.data.remote` to convert trait data from the external TRY database into the tabular format required by the PEcAn meta-analysis module (#3717).
 - Add function `qsub_sda()` for submitting SDA batch jobs by splitting a large number of sites into multiple small groups of sites (#3634).
diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R
index 2b7b8ef0e27..4d39fd3c37e 100644
--- a/base/workflow/R/run.write.configs.R
+++ b/base/workflow/R/run.write.configs.R
@@ -36,9 +36,9 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
   # number of (pft, trait, quantile) combinations
   if (!is.null(input_design) && "ensemble" %in% names(settings)) {
     if (nrow(input_design) != ensemble.size) {
-      stop(
+      PEcAn.logger::logger.error(
         "input_design has ", nrow(input_design), " rows, but ensemble.size is ",
-        ensemble.size, ".The design matrix must have exactly one row per run."
+        ensemble.size, ". The design matrix must have exactly one row per run."
       )
     }
   }
@@ -137,6 +137,26 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
         !is.null(input_design) && 
         "param" %in% colnames(input_design)) {
       trait_sample_indices <- input_design[["param"]]
+      if (is.null(trait.samples) || length(trait.samples) == 0) {
+        PEcAn.logger::logger.error(
+          "samples.Rdata does not contain trait.samples required for input_design$param"
+        )
+      }
+      trait_sample_indices <- as.integer(trait_sample_indices)
+      if (any(is.na(trait_sample_indices)) || any(trait_sample_indices < 1L)) {
+        PEcAn.logger::logger.error("input_design$param must contain positive integer indices")
+      }
+      parameter_bank_size <- PEcAn.uncertainty::trait_sample_bank_size(trait.samples)
+      if (parameter_bank_size == 0L) {
+        PEcAn.logger::logger.error(
+          "samples.Rdata does not contain usable trait.samples required for input_design$param"
+        )
+      }
+      if (any(trait_sample_indices > parameter_bank_size)) {
+        PEcAn.logger::logger.error(
+          "input_design$param includes indices beyond the available parameter sample bank"
+        )
+      }
       ensemble.samples <- list()
       for (pft in names(trait.samples)) {
         pft_traits <- trait.samples[[pft]]
@@ -148,6 +168,10 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
         )
         names(ensemble.samples[[pft]]) <- names(pft_traits)
       }
+    } else if ("ensemble" %in% names(settings) && !is.null(input_design)) {
+      PEcAn.logger::logger.error(
+        "input_design must include a `param` column selecting rows from trait.samples"
+      )
     } else {
       # use pre-generated samples
       ensemble.samples <- existing_data$ensemble.samples
diff --git a/base/workflow/tests/testthat/test.run.write.configs.sobol.R b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
new file mode 100644
index 00000000000..c66b2173a25
--- /dev/null
+++ b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
@@ -0,0 +1,83 @@
+test_that("runModule.run.write.configs uses input_design row count", {
+  settings <- PEcAn.settings::Settings(
+    ensemble = list(
+      size = 3,
+      samplingspace = list(parameters = list(method = "uniform"))
+    ),
+    database = list(bety = list(write = FALSE)),
+    pfts = list(list(posterior.files = "post.distns.Rdata"))
+  )
+  input_design <- data.frame(param = seq_len(5))
+  captured <- new.env(parent = emptyenv())
+
+  mockery::stub(
+    runModule.run.write.configs,
+    "PEcAn.workflow::run.write.configs",
+    function(settings, ensemble.size, input_design, write, posterior.files, overwrite) {
+      captured$ensemble.size <- ensemble.size
+      captured$input_design <- input_design
+      list(
+        ensemble = list(ensemble.id = 123),
+        pfts = settings$pfts
+      )
+    }
+  )
+
+  result <- runModule.run.write.configs(
+    settings = settings,
+    input_design = input_design
+  )
+
+  expect_equal(captured$ensemble.size, 5)
+  expect_identical(captured$input_design, input_design)
+  expect_equal(result$ensemble$ensemble.id, 123)
+})
+
+test_that("run.write.configs validates param against the shortest trait sample bank", {
+  withr::with_tempdir({
+    trait.samples <- list(
+      pftA = list(SLA = seq_len(4)),
+      pftB = list(SLA = seq_len(3))
+    )
+    sa.samples <- list()
+    runs.samples <- list()
+    env.samples <- list()
+    ensemble.samples <- list()
+    save(
+      trait.samples,
+      sa.samples,
+      runs.samples,
+      env.samples,
+      ensemble.samples,
+      file = file.path(getwd(), "samples.Rdata")
+    )
+
+    settings <- list(
+      outdir = getwd(),
+      database = list(),
+      model = list(type = "FAKE"),
+      ensemble = list(size = 2),
+      pfts = list(
+        list(name = "pftA", posteriorid = NULL),
+        list(name = "pftB", posteriorid = NULL)
+      )
+    )
+    run_write_configs <- PEcAn.workflow::run.write.configs
+    mockery::stub(
+      run_write_configs,
+      "PEcAn.logger::logger.error",
+      function(...) stop(paste(...), call. = FALSE)
+    )
+
+    expect_error(
+      run_write_configs(
+        settings = settings,
+        ensemble.size = 2,
+        input_design = data.frame(param = c(1L, 4L)),
+        write = FALSE,
+        overwrite = FALSE
+      ),
+      "input_design\\$param includes indices beyond the available parameter sample bank"
+    )
+  })
+})
diff --git a/modules/uncertainty/DESCRIPTION b/modules/uncertainty/DESCRIPTION
index d9577316963..d7973645550 100644
--- a/modules/uncertainty/DESCRIPTION
+++ b/modules/uncertainty/DESCRIPTION
@@ -43,7 +43,8 @@ Imports:
     purrr,
     randtoolbox,
     rlang,
-    sensitivity
+    sensobol,
+    tibble
 Suggests:
     mockery,
     testthat (>= 1.0.2),
diff --git a/modules/uncertainty/NAMESPACE b/modules/uncertainty/NAMESPACE
index c3ee9bc89f2..5fc8e37759f 100644
--- a/modules/uncertainty/NAMESPACE
+++ b/modules/uncertainty/NAMESPACE
@@ -33,6 +33,7 @@ export(sd.var)
 export(sensitivity.analysis)
 export(sensitivity.filename)
 export(spline.truncate)
+export(trait_sample_bank_size)
 export(write.ensemble.configs)
 export(write.sa.configs)
 importFrom(dplyr,"%>%")
diff --git a/modules/uncertainty/NEWS.md b/modules/uncertainty/NEWS.md
index 5a446773800..4614fc8e83a 100644
--- a/modules/uncertainty/NEWS.md
+++ b/modules/uncertainty/NEWS.md
@@ -1,5 +1,7 @@
 # PEcAn.uncertainty 1.9.0.9000
 
+* Switched Sobol global sensitivity workflows to `sensobol`, including explicit design metadata and saved Sobol index outputs.
+* Treat all inputs (met, initial conditions, events, etc..) as independent Sobol factors; removed parent-child filter that previously confounded events with meteorology.
 * run.ensemble.analysis() now respects `settings$modeloutdir` rather than assuming an `out/` folder inside `settings$outdir` (@Akash-paluvai, #3722).
 * Added `generate_OAT_SA_design()` for creating input design matrices for
 sensitivity analysis. This function ensures non-parameter inputs
diff --git a/modules/uncertainty/R/compute_sobol_indices.R b/modules/uncertainty/R/compute_sobol_indices.R
index d565a771f63..e30e62459c0 100644
--- a/modules/uncertainty/R/compute_sobol_indices.R
+++ b/modules/uncertainty/R/compute_sobol_indices.R
@@ -1,45 +1,145 @@
 #' Compute Sobol indices from a finished PEcAn run
 #'
-#' Loads model outputs from a Sobol ensemble, calculates summary
-#' statistics for a chosen variable, feeds them to \code{sensitivity::tell()},
-#' and returns the updated Sobol object.
+#' Loads standardized ensemble output for a Sobol run, computes first-order and
+#' total-order Sobol indices with \code{sensobol}, and saves both the Sobol design
+#' metadata and the computed indices using PEcAn-style filenames.
 #'
-#' @param outdir     PEcAn run output directory that contains runs.txt
-#' @param sobol_obj  object produced by PEcAn.uncertainty::generate_joint_ensemble_design()
-#' @param var        Variable name to summarise (default "GPP").
-#' @param stat_fun   Summary statistic applied to var default mean .
+#' First-order indices quantify the share of output variance attributable to a
+#' factor alone, while total-order indices summarize the full contribution of
+#' that factor including its interactions with other factors. PEcAn uses
+#' \code{sensobol} for these variance-based estimators; see Saltelli et al. (2008)
+#' for methodological background and Puy et al. (2022) for package details.
 #'
-#' @return           sobol_obj 
-#' .
+#' This function handles one output variable at a time. To compute indices for
+#' multiple variables, call it in a loop (see examples).
+#'
+#' @param outdir PEcAn run output directory containing \code{ensemble.output.*.Rdata}
+#'   files.
+#' @param sobol_obj object produced by
+#'   \code{PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)}.
+#' @param var Variable name to summarize (default \code{"GPP"}).
+#'
+#' @examples
+#' \dontrun{
+#'   # single variable
+#'   result <- compute_sobol_indices(outdir, sobol_obj, var = "GPP")
+#'
+#'   # multiple variables
+#'   vars <- c("GPP", "NPP", "TotSoilCarb")
+#'   all_results <- purrr::map_dfr(vars, function(v) {
+#'     compute_sobol_indices(outdir, sobol_obj, var = v) |>
+#'       dplyr::mutate(variable = v)
+#'   })
+#' }
+#'
+#' @return A tibble of Sobol first-order and total-order indices with attached
+#'   factor metadata.
+#' @references Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni,
+#'   J., Gatelli, D., et al. (2008). Global Sensitivity Analysis: The Primer.
+#'   John Wiley & Sons.
+#'
+#'   Puy, A., Lo Piano, S., Saltelli, A., and Levin, S. A. (2022).
+#'   sensobol: An R Package to Compute Variance-Based Sensitivity Indices.
+#'   Journal of Statistical Software, 102(5), 1-37.
+#'   \doi{10.18637/jss.v102.i05}
 #' @export
 compute_sobol_indices <- function(outdir,
                                   sobol_obj,
-                                  var = "GPP",
-                                  stat_fun = mean) {
- 
- 
- 
-  runs_file <- file.path(outdir, "runs.txt")
-  if (!file.exists(runs_file)) {
-    PEcAn.logger::logger.error("runs.txt not found in ", outdir)
-   }
-  run_ids <- readLines(runs_file) 
-  
-
-  
-  # Load outputs and compute response vector y
-  y <- vapply(run_ids, function(rid) {
-    fpath <- file.path(outdir, rid)
-    out   <- PEcAn.utils::read.output(runid = rid, outdir = fpath)
-    if (!is.list(out) || !var %in% names(out)) {
-      PEcAn.logger::logger.error("Variable '", var, "' missing in output for run ", rid)
-    }
-    stat_fun(out[[var]], na.rm = TRUE)
-  }, numeric(1))
-  
-  # Compute Sobol indices
-  sobol_obj <-sensitivity::tell(sobol_obj, y)
-  
-  # Return the updated object
-  return(invisible(sobol_obj))
+                                  var = "GPP") {
+  if (is.null(sobol_obj$backend) || sobol_obj$backend != "sensobol") {
+    PEcAn.logger::logger.error(
+      "compute_sobol_indices expects a sensobol design object returned by ",
+      "generate_joint_ensemble_design(..., sobol = TRUE)"
+    )
+  }
+
+  output_files <- list.files(
+    outdir,
+    pattern = "^ensemble\\.output\\..*\\.Rdata$",
+    full.names = TRUE
+  )
+  if (length(output_files) == 0) {
+    PEcAn.logger::logger.error("No ensemble.output.*.Rdata files found in ", outdir)
+  }
+
+  output_var <- vapply(
+    strsplit(basename(output_files), "\\."),
+    function(x) if (length(x) >= 4) x[[4]] else NA_character_,
+    character(1)
+  )
+  matched_files <- output_files[output_var == var]
+
+  if (length(matched_files) == 0) {
+    PEcAn.logger::logger.error(
+      "No standardized ensemble output found for variable '", var,
+      "' in ", outdir
+    )
+  }
+  if (length(matched_files) > 1) {
+    PEcAn.logger::logger.error(
+      "Multiple standardized ensemble outputs found for variable '", var,
+      "' in ", outdir, ". Please keep only one matching file."
+    )
+  }
+
+  output_file <- matched_files[[1]]
+  output_env <- new.env(parent = emptyenv())
+  load(output_file, envir = output_env)
+  if (is.null(output_env$ensemble.output)) {
+    PEcAn.logger::logger.error(
+      "Object `ensemble.output` missing from standardized output file ",
+      output_file
+    )
+  }
+
+  y <- as.numeric(unlist(output_env$ensemble.output, use.names = FALSE))
+  expected_length <- sobol_obj$N * (length(sobol_obj$params) + 2L)
+  if (length(y) != expected_length) {
+    PEcAn.logger::logger.error(
+      "Standardized ensemble output has ", length(y),
+      " values but expected ", expected_length,
+      " for Sobol design size N = ", sobol_obj$N
+    )
+  }
+
+  sobol_indices_result <- sensobol::sobol_indices(
+    matrices = sobol_obj$matrices,
+    Y = y,
+    N = sobol_obj$N,
+    params = sobol_obj$params,
+    first = sobol_obj$first,
+    total = sobol_obj$total,
+    order = "first",
+    boot = FALSE
+  )
+
+  sobol_results <- tibble::as_tibble(sobol_indices_result$results)
+  if (!is.null(sobol_obj$factor_metadata)) {
+    factor_metadata <- tibble::as_tibble(sobol_obj$factor_metadata) |>
+      dplyr::rename(parameters = .data$factor)
+    sobol_results <- dplyr::left_join(
+      sobol_results,
+      factor_metadata,
+      by = "parameters"
+    )
+  }
+
+  sobol_design <- sobol_obj
+  save(
+    sobol_design,
+    file = file.path(
+      outdir,
+      sub("^ensemble\\.output\\.", "sobol.design.", basename(output_file))
+    )
+  )
+  save(
+    sobol_indices_result,
+    sobol_results,
+    file = file.path(
+      outdir,
+      sub("^ensemble\\.output\\.", "sobol.indices.", basename(output_file))
+    )
+  )
+
+  return(sobol_results)
 }
diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 0d1a9115ab6..854c5b18dd5 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -598,21 +598,20 @@ input.ens.gen <- function(settings, ensemble_size, input, method = "sampling", p
   if (input == "parameters") return(NULL)
 
   input_path <- settings$run$inputs[[tolower(input)]]$path
-  if (is.null(input_path)) {
-    PEcAn.logger::logger.error(
-      "No paths found for input", sQuote(input), "in settings$run$inputs"
-    )
+  if (is.null(input_path) || length(input_path) == 0) {
+    PEcAn.logger::logger.error("Input ", sQuote(input), " has no paths specified")
   }
 
   if (!is.null(parent_ids)) {
     samples$ids <- parent_ids$ids
-    out.of.sample.size <- length(samples$ids[samples$ids > length(input_path)])
-    # sample for those that our outside the param size -
-    # for example, parent id may send id number 200 but we have only 100 sample for param
-    samples$ids[samples$ids %in% out.of.sample.size] <- sample(
-      seq_along(input_path),
-      out.of.sample.size,
-      replace = TRUE)
+    overflow_positions <- which(samples$ids > length(input_path))
+    if (length(overflow_positions) > 0) {
+      samples$ids[overflow_positions] <- sample(
+        seq_along(input_path),
+        length(overflow_positions),
+        replace = TRUE
+      )
+    }
   } else if (tolower(method) == "sampling") {
     samples$ids <- sample(
       seq_along(input_path),
@@ -623,6 +622,14 @@ input.ens.gen <- function(settings, ensemble_size, input, method = "sampling", p
       seq_along(input_path),
       length.out = ensemble_size)
   }
+  samples$ids <- as.integer(samples$ids)
+  if (any(is.na(samples$ids)) ||
+      any(samples$ids < 1L) ||
+      any(samples$ids > length(input_path))) {
+    PEcAn.logger::logger.error(
+      "Invalid sampled ids generated for input ", sQuote(input)
+    )
+  }
   #using the sample ids
   samples$samples <- input_path[samples$ids]
 
diff --git a/modules/uncertainty/R/generate_joint_ensemble_design.R b/modules/uncertainty/R/generate_joint_ensemble_design.R
index c94b0525127..5967ee1ce8a 100644
--- a/modules/uncertainty/R/generate_joint_ensemble_design.R
+++ b/modules/uncertainty/R/generate_joint_ensemble_design.R
@@ -1,122 +1,196 @@
+#' Minimum trait sample bank size across PFTs
+#'
+#' Returns the smallest number of trait samples available across all PFTs
+#' and traits. Used to verify that the parameter bank is large enough for
+#' a Sobol or ensemble design.
+#'
+#' @param trait.samples named list of PFT trait sample lists, as stored
+#'   in \code{samples.Rdata}.
+#' @return integer, minimum bank size (0 if empty)
+#' @keywords internal
+#' @export
+trait_sample_bank_size <- function(trait.samples) {
+  if (is.null(trait.samples) || length(trait.samples) == 0) {
+    return(0L)
+  }
+
+  bank_sizes <- unlist(
+    purrr::map(trait.samples, function(pft_traits) {
+      if (is.null(pft_traits) || length(pft_traits) == 0) {
+        return(integer(0))
+      }
+      purrr::map_int(pft_traits, function(trait_values) {
+        if (is.null(trait_values) || length(trait_values) == 0) {
+          return(NA_integer_)
+        }
+        as.integer(length(trait_values))
+      })
+    }),
+    use.names = FALSE
+  )
+
+  bank_sizes <- bank_sizes[!is.na(bank_sizes) & bank_sizes > 0L]
+  if (length(bank_sizes) == 0) {
+    return(0L)
+  }
+
+  as.integer(min(bank_sizes))
+}
+
+.sobol_parameter_bank_size <- function(samples_file) {
+  if (!file.exists(samples_file)) {
+    return(0L)
+  }
+
+  samples <- new.env(parent = emptyenv())
+  load(samples_file, envir = samples)
+
+  if (is.null(samples$trait.samples) || length(samples$trait.samples) == 0) {
+    return(0L)
+  }
+
+  trait_sample_bank_size(samples$trait.samples)
+}
+
+.map_sobol_to_indices <- function(x, size) {
+  indices <- floor(stats::qunif(x, min = 1, max = size + 1))
+  as.integer(pmin(indices, size))
+}
+
 #' Generate joint ensemble design for parameter sampling
+#'
 #' Creates a joint ensemble design that maintains parameter correlations across
 #' all sites in a multi-site run. This function generates sample indices that
 #' are shared across sites to ensure consistent parameter sampling.
 #'
-#' @details
-#' Note on internal dependencies
-#'
-#' If samples.Rdata doesn't exist we call get.parameter.samples(), which loads
-#' parameter distributions.
-#'
-#' In practice it:
-#' - uses pft$posterior.files directly when it is defined (an Rdata file with
-#'   post.distns or prior.distns),
-#' - otherwise figures out an output directory from pft$outdir or, if needed,
-#'   via pft$posteriorid in the database,
-#' - then looks in that directory for post.distns.Rdata, falling back to
-#'   prior.distns.Rdata,
-#' - and, for MCMC posteriors, looks up trait.mcmc*.Rdata linked to the same
-#'   posteriorid or a trait.mcmc.Rdata file in that directory.
+#' When \code{sobol = TRUE}, every input listed in
+#' \code{settings$ensemble$samplingspace} (other than \code{parameters})
+#' becomes an independent Sobol factor. This allows variance-based
+#' sensitivity analysis to attribute output variance to each source
+#' (parameters, met, initial conditions, events, etc.) independently.
 #'
-#' Difference from generate_OAT_SA_design: This function samples inputs
-#' randomly or quasi-randomly, while generate_OAT_SA_design holds all
-#' non-parameter inputs constant to isolate parameter effects.
-#'
-#' @param settings PEcAn settings object. This function directly uses:
-#'   \itemize{
-#'     \item \code{settings$outdir} - Output directory path for samples.Rdata
-#'     \item \code{settings$pfts} - List of PFTs (extracts \code{posterior.files})
-#'     \item \code{settings$ensemble$samplingspace} - Input sampling configuration
-#'     \item \code{settings$run$inputs} - Input paths for each input type
-#'   }
-#'   When samples.Rdata doesn't exist, settings is passed to
-#'   \code{\link{get.parameter.samples}} which additionally requires:
-#'   \itemize{
-#'     \item \code{settings$ensemble} - Ensemble configuration
-#'     \item \code{settings$database$bety} - Database connection (optional)
-#'     \item \code{settings$host$name} - Host name for dbfile.check (optional)
-#'   }
+#' @param settings A PEcAn settings object containing ensemble configuration.
 #' @param ensemble_size Integer specifying the number of ensemble members.
-#'   The input_design is generated once for the entire model run. You might
-#'   want to recycle existing ensemble_samples when splitting larger runs
-#'   into smaller jobs while keeping the same parameters.
-#' @param sobol Logical. If TRUE, returns a \code{sensitivity::soboljansen}
-#'   object for Sobol sensitivity analysis.
+#'   When \code{sobol = TRUE}, this is the Sobol base sample size \code{N}, not
+#'   the expanded number of model runs (which will be \code{N * (k + 2)} for
+#'   \code{k} independent factors).
+#' @param sobol Logical, generate a variance-based Sobol design using
+#'   \code{sensobol}.
 #'
-#' @return A list containing ensemble samples and indices.
-#'   If \code{sobol = FALSE}, returns \code{list(X = design_matrix)}.
-#'   If \code{sobol = TRUE}, returns a \code{sensitivity::soboljansen()}
-#'   result object with the design matrix in \code{$X} plus additional
-#'   components for Sobol index calculations.
+#' @return A list with component \code{X}, a data frame design matrix
+#'   describing PEcAn parameter and sampled-input indices. If \code{sobol = TRUE},
+#'   the list also includes the metadata needed by
+#'   \code{\link{compute_sobol_indices}}: \code{N}, \code{params},
+#'   \code{backend}, \code{matrices}, \code{first}, \code{total}, and
+#'   \code{factor_metadata}.
 #'
 #' @export
 
 generate_joint_ensemble_design <- function(settings,
                                            ensemble_size,
                                            sobol = FALSE) {
-  if (sobol) {
-    ensemble_size <- as.numeric(ensemble_size) * 2
-  }
   ens.sample.method <- settings$ensemble$samplingspace$parameters$method
   design_list <- list()
   sampled_inputs <- list()
-  posterior.files <- settings$pfts %>%
+  posterior.files <- settings$pfts |>
     purrr::map_chr("posterior.files", .default = NA_character_)
   samp <- settings$ensemble$samplingspace
-  parents <- lapply(samp, "[[", "parent")
-  order <- names(samp)[
-    lapply(parents, function(tr) which(names(samp) %in% tr)) %>%
-      unlist()
-  ]
-  samp.ordered <- samp[c(order, names(samp)[!(names(samp) %in% order)])]
-
-  # loop over inputs.
-  for (i in seq_along(samp.ordered)) {
-    input_tag <- names(samp.ordered)[i]
-    parent_name <- samp.ordered[[i]]$parent
-
-    parent_ids <- if (!is.null(parent_name)) {
-      sampled_inputs[[parent_name]]
-    } else {
-      NULL
+
+  if (sobol) {
+    # every input in samplingspace (except parameters) is an independent factor
+    input_names <- setdiff(names(samp), "parameters")
+    sobol_factors <- c("param", input_names)
+
+    total_runs <- as.integer(ensemble_size) * (length(sobol_factors) + 2L)
+    samples_file <- file.path(settings$outdir, "samples.Rdata")
+    if (.sobol_parameter_bank_size(samples_file) < total_runs) {
+      PEcAn.uncertainty::get.parameter.samples(
+        settings = settings,
+        ensemble.size = total_runs,
+        posterior.files = posterior.files,
+        ens.sample.method = ens.sample.method
+      )
+    }
+
+    sobol_design <- sensobol::sobol_matrices(
+      matrices = c("A", "B", "AB"),
+      N = as.integer(ensemble_size),
+      params = sobol_factors,
+      order = "first",
+      type = "QRN"
+    )
+    sobol_design <- as.data.frame(sobol_design)
+
+    # map param column to trait bank indices
+    sobol_indices <- list()
+    sobol_indices[["param"]] <- .map_sobol_to_indices(
+      sobol_design[["param"]],
+      total_runs
+    )
+    sampled_inputs[["parameters"]] <- list(ids = sobol_indices[["param"]])
+
+    # map each input to its available paths
+    for (input_tag in input_names) {
+      input_paths <- settings$run$inputs[[tolower(input_tag)]]$path
+      if (is.null(input_paths) || length(input_paths) == 0) {
+        PEcAn.logger::logger.error(
+          "Input ", sQuote(input_tag), " has no paths specified"
+        )
+      }
+      sobol_indices[[input_tag]] <- .map_sobol_to_indices(
+        sobol_design[[input_tag]],
+        length(input_paths)
+      )
+      sampled_inputs[[input_tag]] <- list(ids = sobol_indices[[input_tag]])
+      design_list[[input_tag]] <- sobol_indices[[input_tag]]
     }
 
+    design_list[["param"]] <- sobol_indices[["param"]]
+    design_matrix <- tibble::as_tibble(design_list)
+
+    factor_metadata <- tibble::tibble(
+      factor = sobol_factors,
+      source_type = sobol_factors,
+      source_tag = ifelse(
+        sobol_factors == "param", NA_character_, sobol_factors
+      )
+    )
+
+    return(list(
+      X = design_matrix,
+      N = as.integer(ensemble_size),
+      params = sobol_factors,
+      backend = "sensobol",
+      matrices = c("A", "B", "AB"),
+      first = "saltelli",
+      total = "jansen",
+      factor_metadata = factor_metadata
+    ))
+  }
+
+  # non-Sobol path: simple sequential or sampled design
+  sampled_inputs[["parameters"]] <- list(ids = seq_len(ensemble_size))
+  for (input_tag in setdiff(names(samp), "parameters")) {
     input_result <- PEcAn.uncertainty::input.ens.gen(
       settings = settings,
       ensemble_size = ensemble_size,
       input = input_tag,
-      method = samp.ordered[[i]]$method,
-      parent_ids = parent_ids
+      method = samp[[input_tag]]$method
     )
-
-    sampled_inputs[[input_tag]] <- input_result$ids
+    sampled_inputs[[input_tag]] <- input_result
     design_list[[input_tag]] <- input_result$ids
   }
-  # Sample parameters if we don't have it.
+
   if (!file.exists(file.path(settings$outdir, "samples.Rdata"))) {
     PEcAn.uncertainty::get.parameter.samples(
       settings,
       ensemble.size = ensemble_size,
       posterior.files,
-      ens.sample.method)
+      ens.sample.method
+    )
   }
-  # Here we assumed the length of parameters is identical to the ensemble size.
-  # TODO: detect if they are identical. If not, we will need to resample the 
-  # parameters with replacement.
-  design_list[["param"]] <- seq_len(ensemble_size)
-  design_matrix <- data.frame(design_list)
 
-  if (sobol) {
-    half <- floor(ensemble_size / 2)
-    X1 <- design_matrix[1:half, ]
-    X2 <- design_matrix[(half + 1):ensemble_size, ]
-    sobol_obj <- sensitivity::soboljansen(model = NULL, X1 = X1, X2 = X2)
-    return(sobol_obj)
-  }
-  # This ensures that regardless of whether the sobol or non-sobol version is called 
-  # that the output is a list that includes the design as X. In the sobol version the 
-  # list includes additional info beyond just X that's required by the function that 
-  # does the sobol index calculations, but not required to do the runs themselves.
+  design_list[["param"]] <- sampled_inputs[["parameters"]]$ids
+  design_matrix <- tibble::as_tibble(design_list)
   return(list(X = design_matrix))
-}
\ No newline at end of file
+}
diff --git a/modules/uncertainty/man/compute_sobol_indices.Rd b/modules/uncertainty/man/compute_sobol_indices.Rd
index 4ebf55ece50..b05e6c1f9a0 100644
--- a/modules/uncertainty/man/compute_sobol_indices.Rd
+++ b/modules/uncertainty/man/compute_sobol_indices.Rd
@@ -4,23 +4,57 @@
 \alias{compute_sobol_indices}
 \title{Compute Sobol indices from a finished PEcAn run}
 \usage{
-compute_sobol_indices(outdir, sobol_obj, var = "GPP", stat_fun = mean)
+compute_sobol_indices(outdir, sobol_obj, var = "GPP")
 }
 \arguments{
-\item{outdir}{PEcAn run output directory that contains runs.txt}
+\item{outdir}{PEcAn run output directory containing \code{ensemble.output.*.Rdata}
+files.}
 
-\item{sobol_obj}{object produced by PEcAn.uncertainty::generate_joint_ensemble_design()}
+\item{sobol_obj}{object produced by
+\code{PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)}.}
 
-\item{var}{Variable name to summarise (default "GPP").}
-
-\item{stat_fun}{Summary statistic applied to var default mean .}
+\item{var}{Variable name to summarize (default \code{"GPP"}).}
 }
 \value{
-sobol_obj 
-.
+A tibble of Sobol first-order and total-order indices with attached
+  factor metadata.
 }
 \description{
-Loads model outputs from a Sobol ensemble, calculates summary
-statistics for a chosen variable, feeds them to \code{sensitivity::tell()},
-and returns the updated Sobol object.
+Loads standardized ensemble output for a Sobol run, computes first-order and
+total-order Sobol indices with \code{sensobol}, and saves both the Sobol design
+metadata and the computed indices using PEcAn-style filenames.
+}
+\details{
+First-order indices quantify the share of output variance attributable to a
+factor alone, while total-order indices summarize the full contribution of
+that factor including its interactions with other factors. PEcAn uses
+\code{sensobol} for these variance-based estimators; see Saltelli et al. (2008)
+for methodological background and Puy et al. (2022) for package details.
+
+This function handles one output variable at a time. To compute indices for
+multiple variables, call it in a loop (see examples).
+}
+\examples{
+\dontrun{
+  # single variable
+  result <- compute_sobol_indices(outdir, sobol_obj, var = "GPP")
+
+  # multiple variables
+  vars <- c("GPP", "NPP", "TotSoilCarb")
+  all_results <- purrr::map_dfr(vars, function(v) {
+    compute_sobol_indices(outdir, sobol_obj, var = v) |>
+      dplyr::mutate(variable = v)
+  })
+}
+
+}
+\references{
+Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni,
+  J., Gatelli, D., et al. (2008). Global Sensitivity Analysis: The Primer.
+  John Wiley & Sons.
+
+  Puy, A., Lo Piano, S., Saltelli, A., and Levin, S. A. (2022).
+  sensobol: An R Package to Compute Variance-Based Sensitivity Indices.
+  Journal of Statistical Software, 102(5), 1-37.
+  \doi{10.18637/jss.v102.i05}
 }
diff --git a/modules/uncertainty/man/generate_joint_ensemble_design.Rd b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
index 77ccc765227..a1e3406bee7 100644
--- a/modules/uncertainty/man/generate_joint_ensemble_design.Rd
+++ b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
@@ -2,67 +2,38 @@
 % Please edit documentation in R/generate_joint_ensemble_design.R
 \name{generate_joint_ensemble_design}
 \alias{generate_joint_ensemble_design}
-\title{Generate joint ensemble design for parameter sampling
-Creates a joint ensemble design that maintains parameter correlations across
-all sites in a multi-site run. This function generates sample indices that
-are shared across sites to ensure consistent parameter sampling.}
+\title{Generate joint ensemble design for parameter sampling}
 \usage{
 generate_joint_ensemble_design(settings, ensemble_size, sobol = FALSE)
 }
 \arguments{
-\item{settings}{PEcAn settings object. This function directly uses:
-\itemize{
-  \item \code{settings$outdir} - Output directory path for samples.Rdata
-  \item \code{settings$pfts} - List of PFTs (extracts \code{posterior.files})
-  \item \code{settings$ensemble$samplingspace} - Input sampling configuration
-  \item \code{settings$run$inputs} - Input paths for each input type
-}
-When samples.Rdata doesn't exist, settings is passed to
-\code{\link{get.parameter.samples}} which additionally requires:
-\itemize{
-  \item \code{settings$ensemble} - Ensemble configuration
-  \item \code{settings$database$bety} - Database connection (optional)
-  \item \code{settings$host$name} - Host name for dbfile.check (optional)
-}}
+\item{settings}{A PEcAn settings object containing ensemble configuration.}
 
 \item{ensemble_size}{Integer specifying the number of ensemble members.
-The input_design is generated once for the entire model run. You might
-want to recycle existing ensemble_samples when splitting larger runs
-into smaller jobs while keeping the same parameters.}
+When \code{sobol = TRUE}, this is the Sobol base sample size \code{N}, not
+the expanded number of model runs (which will be \code{N * (k + 2)} for
+\code{k} independent factors).}
 
-\item{sobol}{Logical. If TRUE, returns a \code{sensitivity::soboljansen}
-object for Sobol sensitivity analysis.}
+\item{sobol}{Logical, generate a variance-based Sobol design using
+\code{sensobol}.}
 }
 \value{
-A list containing ensemble samples and indices.
-  If \code{sobol = FALSE}, returns \code{list(X = design_matrix)}.
-  If \code{sobol = TRUE}, returns a \code{sensitivity::soboljansen()}
-  result object with the design matrix in \code{$X} plus additional
-  components for Sobol index calculations.
+A list with component \code{X}, a data frame design matrix
+  describing PEcAn parameter and sampled-input indices. If \code{sobol = TRUE},
+  the list also includes the metadata needed by
+  \code{\link{compute_sobol_indices}}: \code{N}, \code{params},
+  \code{backend}, \code{matrices}, \code{first}, \code{total}, and
+  \code{factor_metadata}.
 }
 \description{
-Generate joint ensemble design for parameter sampling
 Creates a joint ensemble design that maintains parameter correlations across
 all sites in a multi-site run. This function generates sample indices that
 are shared across sites to ensure consistent parameter sampling.
 }
 \details{
-Note on internal dependencies
-
-If samples.Rdata doesn't exist we call get.parameter.samples(), which loads
-parameter distributions.
-
-In practice it:
-- uses pft$posterior.files directly when it is defined (an Rdata file with
-  post.distns or prior.distns),
-- otherwise figures out an output directory from pft$outdir or, if needed,
-  via pft$posteriorid in the database,
-- then looks in that directory for post.distns.Rdata, falling back to
-  prior.distns.Rdata,
-- and, for MCMC posteriors, looks up trait.mcmc*.Rdata linked to the same
-  posteriorid or a trait.mcmc.Rdata file in that directory.
-
-Difference from generate_OAT_SA_design: This function samples inputs
-randomly or quasi-randomly, while generate_OAT_SA_design holds all
-non-parameter inputs constant to isolate parameter effects.
+When \code{sobol = TRUE}, every input listed in
+\code{settings$ensemble$samplingspace} (other than \code{parameters})
+becomes an independent Sobol factor. This allows variance-based
+sensitivity analysis to attribute output variance to each source
+(parameters, met, initial conditions, events, etc.) independently.
 }
diff --git a/modules/uncertainty/man/trait_sample_bank_size.Rd b/modules/uncertainty/man/trait_sample_bank_size.Rd
new file mode 100644
index 00000000000..9101f8988d6
--- /dev/null
+++ b/modules/uncertainty/man/trait_sample_bank_size.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/generate_joint_ensemble_design.R
+\name{trait_sample_bank_size}
+\alias{trait_sample_bank_size}
+\title{Minimum trait sample bank size across PFTs}
+\usage{
+trait_sample_bank_size(trait.samples)
+}
+\arguments{
+\item{trait.samples}{named list of PFT trait sample lists, as stored
+in \code{samples.Rdata}.}
+}
+\value{
+integer, minimum bank size (0 if empty)
+}
+\description{
+Returns the smallest number of trait samples available across all PFTs
+and traits. Used to verify that the parameter bank is large enough for
+a Sobol or ensemble design.
+}
+\keyword{internal}
diff --git a/modules/uncertainty/tests/testthat/test.sobol.R b/modules/uncertainty/tests/testthat/test.sobol.R
new file mode 100644
index 00000000000..4ffcea7aa42
--- /dev/null
+++ b/modules/uncertainty/tests/testthat/test.sobol.R
@@ -0,0 +1,217 @@
+make_sobol_settings <- function(outdir) {
+  samples_file <- file.path(outdir, "samples.Rdata")
+  trait.samples <- list(
+    temperate = list(
+      SLA = seq_len(40)
+    )
+  )
+  save(trait.samples, file = samples_file)
+
+  PEcAn.settings::Settings(
+    outdir = outdir,
+    pfts = list(list(name = "temperate", posterior.files = "post.distns.Rdata")),
+    run = list(inputs = list(
+      met = list(path = c("met1", "met2", "met3")),
+      poolinitcond = list(path = c("ic1", "ic2")),
+      events = list(path = c("evt1", "evt2", "evt3"))
+    )),
+    ensemble = list(
+      samplingspace = list(
+        parameters = list(method = "uniform"),
+        met = list(method = "sampling"),
+        poolinitcond = list(method = "looping"),
+        events = list(method = "sampling")
+      )
+    )
+  )
+}
+
+make_mixed_bank_sobol_settings <- function(outdir) {
+  samples_file <- file.path(outdir, "samples.Rdata")
+  trait.samples <- list(
+    hardwood = list(
+      SLA = seq_len(20)
+    ),
+    conifer = list(
+      SLA = seq_len(10)
+    )
+  )
+  save(trait.samples, file = samples_file)
+
+  PEcAn.settings::Settings(
+    outdir = outdir,
+    pfts = list(
+      list(name = "hardwood", posterior.files = "post1.distns.Rdata"),
+      list(name = "conifer", posterior.files = "post2.distns.Rdata")
+    ),
+    run = list(inputs = list()),
+    ensemble = list(
+      samplingspace = list(
+        parameters = list(method = "uniform")
+      )
+    )
+  )
+}
+
+test_that("Sobol design treats all inputs as independent factors", {
+  withr::with_tempdir({
+    settings <- make_sobol_settings(getwd())
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 4,
+      sobol = TRUE
+    )
+
+    # 4 independent factors: param, met, poolinitcond, events
+    # total runs = N * (k + 2) = 4 * (4 + 2) = 24
+    expect_equal(nrow(result$X), 24)
+    expect_equal(result$N, 4)
+    expect_identical(
+      result$params,
+      c("param", "met", "poolinitcond", "events")
+    )
+    expect_identical(result$backend, "sensobol")
+    expect_identical(result$matrices, c("A", "B", "AB"))
+    expect_identical(result$first, "saltelli")
+    expect_identical(result$total, "jansen")
+
+    # all factor columns present
+    expect_true(all(
+      c("param", "met", "poolinitcond", "events") %in% names(result$X)
+    ))
+
+    # events should have independent indices -- not identical to met
+    # (quasi-random design makes exact equality extremely unlikely)
+    expect_false(identical(result$X$events, result$X$met))
+
+    # parameter indices stay within bank range
+    expect_true(all(result$X$param >= 1))
+    expect_true(all(result$X$param <= 24))
+
+    # input indices stay within available paths
+    expect_true(all(result$X$met >= 1 & result$X$met <= 3))
+    expect_true(all(result$X$poolinitcond >= 1 & result$X$poolinitcond <= 2))
+    expect_true(all(result$X$events >= 1 & result$X$events <= 3))
+
+    # factor metadata covers all independent factors
+    expect_true(all(
+      c("factor", "source_type", "source_tag") %in%
+        names(result$factor_metadata)
+    ))
+    expect_identical(
+      result$factor_metadata$source_type,
+      c("param", "met", "poolinitcond", "events")
+    )
+  })
+})
+
+test_that("Non-Sobol design generation remains row-for-row", {
+  withr::with_tempdir({
+    settings <- make_sobol_settings(getwd())
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 5,
+      sobol = FALSE
+    )
+
+    expect_named(result, "X")
+    expect_equal(nrow(result$X), 5)
+    expect_false("backend" %in% names(result))
+  })
+})
+
+test_that("Sobol regenerates parameter bank when any PFT bank is too short", {
+  withr::with_tempdir({
+    settings <- make_mixed_bank_sobol_settings(getwd())
+    captured <- new.env(parent = emptyenv())
+
+    mockery::stub(
+      generate_joint_ensemble_design,
+      "PEcAn.uncertainty::get.parameter.samples",
+      function(settings, ensemble.size, posterior.files, ens.sample.method) {
+        captured$ensemble.size <- ensemble.size
+        captured$posterior.files <- posterior.files
+        captured$ens.sample.method <- ens.sample.method
+        NULL
+      }
+    )
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 5,
+      sobol = TRUE
+    )
+
+    # only param factor here (no inputs in samplingspace)
+    # total = N * (k + 2) = 5 * (1 + 2) = 15
+    expect_equal(captured$ensemble.size, 15)
+    expect_identical(
+      captured$posterior.files,
+      c("post1.distns.Rdata", "post2.distns.Rdata")
+    )
+    expect_identical(captured$ens.sample.method, "uniform")
+    expect_equal(nrow(result$X), 15)
+    expect_true(all(result$X$param >= 1))
+    expect_true(all(result$X$param <= 15))
+  })
+})
+
+test_that("compute_sobol_indices matches direct sensobol results", {
+  withr::with_tempdir({
+    sobol_obj <- list(
+      N = 4L,
+      params = c("param", "met"),
+      backend = "sensobol",
+      matrices = c("A", "B", "AB"),
+      first = "saltelli",
+      total = "jansen",
+      factor_metadata = data.frame(
+        factor = c("param", "met"),
+        source_type = c("param", "met"),
+        source_tag = c(NA_character_, "met"),
+        stringsAsFactors = FALSE
+      )
+    )
+
+    y <- seq_len(16)
+    ensemble.output <- as.list(y)
+    save(ensemble.output, file = file.path(
+      getwd(),
+      "ensemble.output.NOENSEMBLEID.GPP.NA.NA.Rdata"
+    ))
+
+    result <- compute_sobol_indices(
+      outdir = getwd(),
+      sobol_obj = sobol_obj,
+      var = "GPP"
+    )
+
+    expected <- tibble::as_tibble(
+      sensobol::sobol_indices(
+        matrices = c("A", "B", "AB"),
+        Y = y,
+        N = 4L,
+        params = c("param", "met"),
+        first = "saltelli",
+        total = "jansen",
+        order = "first",
+        boot = FALSE
+      )$results
+    )
+
+    expect_equal(
+      result[, c("parameters", "sensitivity", "original")],
+      expected[, c("parameters", "sensitivity", "original")]
+    )
+    expect_true(file.exists(file.path(
+      getwd(),
+      "sobol.design.NOENSEMBLEID.GPP.NA.NA.Rdata"
+    )))
+    expect_true(file.exists(file.path(
+      getwd(),
+      "sobol.indices.NOENSEMBLEID.GPP.NA.NA.Rdata"
+    )))
+  })
+})
diff --git a/modules/uncertainty/tests/testthat/test_ensemble.R b/modules/uncertainty/tests/testthat/test_ensemble.R
index 3d86710807d..81b6678204e 100644
--- a/modules/uncertainty/tests/testthat/test_ensemble.R
+++ b/modules/uncertainty/tests/testthat/test_ensemble.R
@@ -1,5 +1,4 @@
 context("input validation for write.ensemble.configs")
-library(testthat)
 
 # Mock a model write.configs function to avoid model-specific errors
 write.configs.SIPNET <- function(...) TRUE