From d1a2407495d0347762832b12f31fa058b77b9adc Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 4 Mar 2026 23:47:38 -0700
Subject: [PATCH 01/25] Lean SIPNET postprocessing fast-path and runtime
 controls

---
 models/sipnet/DESCRIPTION              |   1 +
 models/sipnet/R/model2netcdf.SIPNET.R  | 568 +++++++++++++++++--------
 models/sipnet/R/write.configs.SIPNET.R | 323 +++++++++-----
 models/sipnet/README.md                |  12 +
 4 files changed, 607 insertions(+), 297 deletions(-)

diff --git a/models/sipnet/DESCRIPTION b/models/sipnet/DESCRIPTION
index 955e2a7403a..21223dca9b1 100644
--- a/models/sipnet/DESCRIPTION
+++ b/models/sipnet/DESCRIPTION
@@ -14,6 +14,7 @@ URL: https://pecanproject.github.io
 BugReports: https://github.com/PecanProject/pecan/issues
 Depends: R (>= 4.1.0)
 Imports:
+    data.table,
     dplyr,
     jsonlite,
     lubridate (>= 1.6.0),
diff --git a/models/sipnet/R/model2netcdf.SIPNET.R b/models/sipnet/R/model2netcdf.SIPNET.R
index 5d7324ef6dc..85c68523f57 100644
--- a/models/sipnet/R/model2netcdf.SIPNET.R
+++ b/models/sipnet/R/model2netcdf.SIPNET.R
@@ -1,15 +1,15 @@
 #' Merge multiple NetCDF files into one
-#' 
+#'
 #' @param files \code{character}. List of filepaths, which should lead to NetCDF files.
 #' @param outfile \code{character}. Output filename of the merged data.
 #' @return A NetCDF file containing all of the merged data.
 #' @examples
 #' \dontrun{
-#' files <- list.files(paste0(system.file(package="processNC"), "/extdata"), 
+#' files <- list.files(paste0(system.file(package="processNC"), "/extdata"),
 #'                     pattern="tas.*\\.nc", full.names=TRUE)
 #' temp <- tempfile(fileext=".nc")
 #' mergeNC(files=files, outfile=temp)
-#' terra::rast(temp) 
+#' terra::rast(temp)
 #' }
 #' @export mergeNC
 #' @name mergeNC
@@ -17,27 +17,82 @@
 mergeNC <- function(
     ##title<< Aggregate data in netCDF files
   files ##<< character vector: names of the files to merge
-  , outfile ##<< character: path to save the results files to. 
+  , outfile ##<< character: path to save the results files to.
 )
   ##description<<
   ## This function aggregates time periods in netCDF files. Basically it is just a
   ## wrapper around the respective cdo function.
 {
-  ##test input
-  #if (system("cdo -V")==0)
-  #  stop('cdo not found. Please install it.')
-  
   ## supply cdo command
-  cdoCmd <- paste('cdo -cat', paste(files, collapse=" "), outfile, sep=' ')
-  
-  ##run command
+  cdoCmd <- paste("cdo -cat", paste(files, collapse = " "), outfile, sep = " ")
+
+  ## run command
   system(cdoCmd)
-  cat(paste('Created file ', outfile, '.\n', sep = ''))
-  
-  ## character string: name of the file created. 
+  cat(paste("Created file ", outfile, ".\n", sep = ""))
+
+  ## character string: name of the file created.
   invisible(outfile)
 }
 
+.sipnet_read_output_table <- function(path, required_cols, optional_cols) {
+  read_header <- function(skip_lines) {
+    try(
+      data.table::fread(
+        path,
+        skip = skip_lines,
+        nrows = 0,
+        data.table = FALSE,
+        check.names = FALSE,
+        showProgress = FALSE
+      ),
+      silent = TRUE
+    )
+  }
+
+  header <- read_header(0L)
+  skip_used <- 0L
+  if (inherits(header, "try-error") || !("year" %in% names(header))) {
+    header <- read_header(1L)
+    skip_used <- 1L
+  }
+
+  if (inherits(header, "try-error") || !("year" %in% names(header))) {
+    PEcAn.logger::logger.error("Unable to parse SIPNET output file header", path)
+  }
+
+  available_cols <- names(header)
+  missing_required <- setdiff(required_cols, available_cols)
+  if (length(missing_required) > 0) {
+    PEcAn.logger::logger.error(
+      "Missing required SIPNET output columns:",
+      paste(missing_required, collapse = ", ")
+    )
+  }
+
+  read_cols <- c(required_cols, intersect(optional_cols, available_cols))
+  int_cols <- intersect(c("year", "day"), read_cols)
+  num_cols <- setdiff(read_cols, int_cols)
+
+  sipnet_output <- try(
+    data.table::fread(
+      path,
+      skip = skip_used,
+      select = read_cols,
+      colClasses = list(integer = int_cols, numeric = num_cols),
+      data.table = FALSE,
+      check.names = FALSE,
+      showProgress = FALSE
+    ),
+    silent = TRUE
+  )
+
+  if (inherits(sipnet_output, "try-error")) {
+    PEcAn.logger::logger.error("Unable to parse SIPNET output file", path)
+  }
+
+  sipnet_output
+}
+
 #--------------------------------------------------------------------------------------------------#
 ##' Convert SIPNET output to netCDF
 ##'
@@ -58,153 +113,230 @@ mergeNC <- function(
 ##' @author Shawn Serbin, Michael Dietze
 model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date, delete.raw = FALSE, revision, prefix = "sipnet.out",
                                 overwrite = FALSE, conflict = FALSE) {
-  ### Read in model output in SIPNET format
+  profile_enabled <- tolower(trimws(Sys.getenv("PECAN_SIPNET_PROFILE", ""))) %in% c("1", "true", "yes", "on")
+  timings <- list()
+  tic <- function() proc.time()[["elapsed"]]
+  toc <- function(stage, started_at) {
+    if (!profile_enabled) {
+      return(invisible(NULL))
+    }
+    timings[[length(timings) + 1]] <<- data.frame(
+      function_name = "model2netcdf.SIPNET",
+      stage = stage,
+      elapsed_seconds = unname(proc.time()[["elapsed"]] - started_at),
+      stringsAsFactors = FALSE
+    )
+    invisible(NULL)
+  }
+
+  write_profile_csv <- trimws(Sys.getenv("PECAN_SIPNET_PROFILE_CSV", ""))
+  write_profile_csv <- if (profile_enabled) {
+    if (tolower(write_profile_csv) %in% c("1", "true", "yes", "on")) {
+      file.path(outdir, "sipnet_model2netcdf.profile.csv")
+    } else if (nzchar(write_profile_csv)) {
+      write_profile_csv
+    } else {
+      ""
+    }
+  } else {
+    ""
+  }
+
+  flush_timings <- function() {
+    if (!profile_enabled || length(timings) == 0) {
+      return(invisible(NULL))
+    }
+    timing_df <- do.call(rbind, timings)
+    summary_txt <- apply(timing_df, 1, function(x) {
+      paste0(x[["function_name"]], " ", x[["stage"]], ": ", signif(as.numeric(x[["elapsed_seconds"]]), 4), " s")
+    })
+    PEcAn.logger::logger.info("SIPNET profiling summary:\n", paste(summary_txt, collapse = "\n"), wrap = FALSE)
+    if (nzchar(write_profile_csv)) {
+      utils::write.table(
+        timing_df,
+        file = write_profile_csv,
+        sep = ",",
+        row.names = FALSE,
+        col.names = !file.exists(write_profile_csv),
+        quote = TRUE,
+        append = file.exists(write_profile_csv)
+      )
+    }
+    invisible(NULL)
+  }
+  on.exit(flush_timings(), add = TRUE)
+
+  t_total <- tic()
   sipnet_out_file <- file.path(outdir, prefix)
-  sipnet_output <- utils::read.table(sipnet_out_file, header = T, skip = 1, sep = "")
-  #sipnet_output_dims <- dim(sipnet_output)
-  
-  ### Determine number of years and output timestep
-  #start.day <- sipnet_output$day[1]
-  num_years <- length(unique(sipnet_output$year))
-  simulation_years <- unique(sipnet_output$year)
-  
-  # get all years that we want data from
+
+  required_cols <- c(
+    "year", "day", "time", "gpp", "rAboveground", "rRoot", "rtot", "rSoil", "nee",
+    "plantWoodC", "plantLeafC", "coarseRootC", "fineRootC", "soil", "litter",
+    "fluxestranspiration", "soilWater", "soilWetnessFrac", "snow"
+  )
+  optional_cols <- c("litterWater", "evapotranspiration", "npp", "woodCreation")
+
+  t_read <- tic()
+  sipnet_output <- .sipnet_read_output_table(sipnet_out_file, required_cols, optional_cols)
+  toc("read_sipnet_output", t_read)
+
+  t_prepare <- tic()
+  simulation_years <- sort(unique(sipnet_output$year))
   year_seq <- seq(lubridate::year(start_date), lubridate::year(end_date))
-  
-  # check that specified years and output years match
+
   if (!all(year_seq %in% simulation_years)) {
-    PEcAn.logger::logger.severe("Years selected for model run and SIPNET output years do not match ")
+    PEcAn.logger::logger.severe("Years selected for model run and SIPNET output years do not match")
+  }
+
+  year_index <- split(seq_len(nrow(sipnet_output)), sipnet_output$year)
+
+  first_year <- simulation_years[1]
+  first_day <- unique(sipnet_output$day[sipnet_output$year == first_year])[1]
+  out_day <- sum(sipnet_output$year == first_year & sipnet_output$day == first_day, na.rm = TRUE)
+  if (!is.finite(out_day) || out_day <= 0) {
+    PEcAn.logger::logger.error("Unable to infer SIPNET output timestep from parsed output")
   }
-  
-  # get number of model timesteps per day
-  # outday is the number of time steps in a day - for example 6 hours would have out_day of 4
-  
-  out_day <- sum(
-    sipnet_output$year == simulation_years[1] &
-      sipnet_output$day == unique(sipnet_output$day)[1],
-    na.rm = TRUE
-  ) # switched to day 2 in case first day is partial
-  
-  
+
   timestep.s <- 86400 / out_day
-  
-  
-  ### Loop over years in SIPNET output to create separate netCDF outputs
-  for (y in year_seq) {
-    #initialize the conflicted as FALSE
-    conflicted <- FALSE
-    conflict <- TRUE    #conflict is set to TRUE to enable the rename of yearly nc file for merging SDA results with sub-annual data
-    #if we have conflicts on this file.
-    if (file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & overwrite == FALSE & conflict == FALSE) {
-      next
-    }else if(file.exists(file.path(outdir, paste(y, "nc", sep = "."))) & conflict){
-      conflicted <- TRUE
-      file.rename(file.path(outdir, paste(y, "nc", sep = ".")), file.path(outdir, "previous.nc"))
+
+  run_dir <- if (grepl("/out/", outdir, fixed = TRUE)) {
+    sub("/out/", "/run/", outdir, fixed = TRUE)
+  } else {
+    sub("/out/?$", "/run/", outdir)
+  }
+  if (identical(run_dir, outdir)) {
+    PEcAn.logger::logger.error("Cannot infer run directory from outdir; expected '/out/' in path", outdir)
+  }
+
+  run_param_file <- file.path(run_dir, "sipnet.param")
+  if (!file.exists(run_param_file)) {
+    PEcAn.logger::logger.error("Missing SIPNET parameter file", run_param_file)
+  }
+
+  param <- utils::read.table(run_param_file, stringsAsFactors = FALSE)
+  id <- which(param[, 1] == "leafCSpWt")
+  SLA <- 1000 / param[id, 2]
+
+  worker_env <- trimws(Sys.getenv("PECAN_SIPNET_NC_WORKERS", ""))
+  if (nzchar(worker_env)) {
+    worker_setting <- suppressWarnings(as.integer(worker_env))
+    if (is.na(worker_setting) || worker_setting < 1L) {
+      PEcAn.logger::logger.error("PECAN_SIPNET_NC_WORKERS must be an integer >= 1")
     }
-    print(paste("---- Processing year: ", y))  # turn on for debugging
-    
-    ## Subset data for processing
-    sub.sipnet.output <- subset(sipnet_output, sipnet_output$year == y)
-    
-    raw_time <- sub.sipnet.output[["time"]] # decimal hours (eg 13.75 = 1:45 PM)
-    doy <- sub.sipnet.output[["day"]] # day of year, not of month
-    hr <- floor(raw_time)
-    minsec <- PEcAn.utils::ud_convert(raw_time - hr, "hour", "min")
-    min <- floor(minsec)
-    sec <- PEcAn.utils::ud_convert(minsec - min, "minute", "second")
-    sub_dates <- strptime(
-      paste(y, doy, hr, min, sec),
-      "%Y %j %H %M %S",
-      tz = "UTC"
-    )
-    sub_dates_cf <- PEcAn.utils::datetime2cf(
-      sub_dates,
-      paste0("days since ", y, "-01-01"),
-      tz = "UTC"
+  } else {
+    detected <- suppressWarnings(as.integer(parallel::detectCores()))
+    if (is.na(detected) || detected < 1L) {
+      detected <- 1L
+    }
+    worker_setting <- max(1L, detected - 1L)
+  }
+
+  if (conflict && !nzchar(Sys.which("cdo"))) {
+    PEcAn.logger::logger.error("'conflict=TRUE' requires cdo in PATH")
+  }
+
+  if (!("litterWater" %in% names(sipnet_output))) {
+    sipnet_output$litterWater <- NA_real_
+  }
+  if (!("woodCreation" %in% names(sipnet_output))) {
+    sipnet_output$woodCreation <- NA_real_
+  }
+
+  if (revision == "unk") {
+    if (!("npp" %in% names(sipnet_output))) {
+      PEcAn.logger::logger.error("SIPNET output missing 'npp' required for revision='unk'")
+    }
+    qle_source <- sipnet_output$npp
+  } else if ("evapotranspiration" %in% names(sipnet_output)) {
+    qle_source <- sipnet_output$evapotranspiration
+  } else if ("npp" %in% names(sipnet_output)) {
+    qle_source <- sipnet_output$npp
+  } else {
+    PEcAn.logger::logger.error("SIPNET output missing both 'evapotranspiration' and 'npp' required for Qle")
+  }
+
+  timestamps_utc <- as.POSIXct(
+    strptime(sprintf("%04d %03d", sipnet_output$year, sipnet_output$day), "%Y %j", tz = "UTC")
+  ) + sipnet_output$time * 3600
+
+  year_origin <- as.POSIXct(
+    paste0(year_seq, "-01-01 00:00:00"),
+    tz = "UTC"
+  )
+  names(year_origin) <- as.character(year_seq)
+
+  output_all <- list(
+    "GPP" = (sipnet_output$gpp * 0.001) / timestep.s,
+    "NPP" = (sipnet_output$gpp * 0.001) / timestep.s -
+      ((sipnet_output$rAboveground * 0.001) / timestep.s + (sipnet_output$rRoot * 0.001) / timestep.s),
+    "TotalResp" = (sipnet_output$rtot * 0.001) / timestep.s,
+    "AutoResp" = (sipnet_output$rAboveground * 0.001) / timestep.s + (sipnet_output$rRoot * 0.001) / timestep.s,
+    "HeteroResp" = ((sipnet_output$rSoil - sipnet_output$rRoot) * 0.001) / timestep.s,
+    "SoilResp" = (sipnet_output$rSoil * 0.001) / timestep.s,
+    "NEE" = (sipnet_output$nee * 0.001) / timestep.s,
+    "AbvGrndWood" = (sipnet_output$plantWoodC * 0.001),
+    "leaf_carbon_content" = (sipnet_output$plantLeafC * 0.001),
+    "TotLivBiom" = (sipnet_output$plantWoodC * 0.001) + (sipnet_output$plantLeafC * 0.001) +
+      (sipnet_output$coarseRootC + sipnet_output$fineRootC) * 0.001,
+    "TotSoilCarb" = (sipnet_output$soil * 0.001) + (sipnet_output$litter * 0.001),
+    "Qle" = (qle_source * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s,
+    "Transp" = (sipnet_output$fluxestranspiration * 10) / timestep.s,
+    "SoilMoist" = (sipnet_output$soilWater * 10),
+    "SoilMoistFrac" = sipnet_output$soilWetnessFrac,
+    "SWE" = (sipnet_output$snow * 10),
+    "litter_carbon_content" = sipnet_output$litter * 0.001,
+    "litter_mass_content_of_water" = sipnet_output$litterWater * 10,
+    "LAI" = (sipnet_output$plantLeafC * 0.001) * SLA,
+    "fine_root_carbon_content" = sipnet_output$fineRootC * 0.001,
+    "coarse_root_carbon_content" = sipnet_output$coarseRootC * 0.001,
+    "GWBI" = (sipnet_output$woodCreation * 0.001) / 86400,
+    "AGB" = (sipnet_output$plantWoodC + sipnet_output$plantLeafC) * 0.001
+  )
+  toc("prepare_conversion_inputs", t_prepare)
+
+  write_year_file <- function(y, row_ids, target_file) {
+    if (length(row_ids) == 0) {
+      return(invisible(NULL))
+    }
+
+    sub_dates_cf <- as.numeric(
+      difftime(timestamps_utc[row_ids], year_origin[[as.character(y)]], units = "days")
     )
-    
-    sub.sipnet.output.dims <- dim(sub.sipnet.output)
     dayfrac <- 1 / out_day
-    
-    # create netCDF time.bounds variable
-    bounds <- array(data=NA, dim=c(length(sub_dates_cf),2))
-    bounds[,1] <- sub_dates_cf
-    bounds[,2] <- bounds[,1]+dayfrac
-    # create time bounds for each timestep in t, t+1; t+1, t+2... format
-    bounds <- round(bounds,4) 
-    
-    ## Setup outputs for netCDF file in appropriate units
-    output       <- list(
-      "GPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s,  # GPP in kgC/m2/s
-      "NPP" = (sub.sipnet.output$gpp * 0.001) / timestep.s - ((sub.sipnet.output$rAboveground *
-                                                                 0.001) / timestep.s + (sub.sipnet.output$rRoot * 0.001) / timestep.s), # NPP in kgC/m2/s. Post SIPNET calculation
-      "TotalResp" = (sub.sipnet.output$rtot * 0.001) / timestep.s,  # Total Respiration in kgC/m2/s
-      "AutoResp" = (sub.sipnet.output$rAboveground * 0.001) / timestep.s + (sub.sipnet.output$rRoot *
-                                                                              0.001) / timestep.s,  # Autotrophic Respiration in kgC/m2/s
-      "HeteroResp" = ((sub.sipnet.output$rSoil - sub.sipnet.output$rRoot) * 0.001) / timestep.s,  # Heterotrophic Respiration in kgC/m2/s
-      "SoilResp" = (sub.sipnet.output$rSoil * 0.001) / timestep.s,  # Soil Respiration in kgC/m2/s
-      "NEE" = (sub.sipnet.output$nee * 0.001) / timestep.s,  # NEE in kgC/m2/s
-      "AbvGrndWood" = (sub.sipnet.output$plantWoodC * 0.001),  # Above ground wood kgC/m2
-      "leaf_carbon_content" = (sub.sipnet.output$plantLeafC * 0.001),  # Leaf C kgC/m2
-      "TotLivBiom" = (sub.sipnet.output$plantWoodC * 0.001) + (sub.sipnet.output$plantLeafC * 0.001) + 
-        (sub.sipnet.output$coarseRootC + sub.sipnet.output$fineRootC) * 0.001, # Total living C kgC/m2
-      "TotSoilCarb" = (sub.sipnet.output$soil * 0.001) + (sub.sipnet.output$litter * 0.001)  # Total soil C kgC/m2
-    )
-    if (revision == "unk") {
-      ## *** NOTE : npp in the sipnet output file is actually evapotranspiration, this is due to a bug in sipnet.c : ***
-      ## *** it says "npp" in the header (written by L774) but the values being written are trackers.evapotranspiration (L806) ***
-      ## evapotranspiration in SIPNET is cm^3 water per cm^2 of area, to convert it to latent heat units W/m2 multiply with :
-      ## 0.01 (cm2m) * 1000 (water density, kg m-3) * latent heat of vaporization (J kg-1)
-      ## latent heat of vaporization is not constant and it varies slightly with temperature, get.lv() returns 2.5e6 J kg-1 by default
-      output[["Qle"]] <- (sub.sipnet.output$npp * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s  # Qle W/m2
-    } else {
-      output[["Qle"]] <- (sub.sipnet.output$evapotranspiration * 10 * PEcAn.data.atmosphere::get.lv()) / timestep.s  # Qle W/m2
-    }
-    output[["Transp"]] <- (sub.sipnet.output$fluxestranspiration * 10) / timestep.s  # Transpiration kgW/m2/s
-    output[["SoilMoist"]] <- (sub.sipnet.output$soilWater * 10)  # Soil moisture kgW/m2
-    output[["SoilMoistFrac"]] <- (sub.sipnet.output$soilWetnessFrac)  # Fractional soil wetness
-    output[["SWE"]] <- (sub.sipnet.output$snow * 10)  # SWE
-    output[["litter_carbon_content"]] <- sub.sipnet.output$litter * 0.001  ## litter kgC/m2
-    output[["litter_mass_content_of_water"]] <- (sub.sipnet.output$litterWater * 10) # Litter water kgW/m2
-    #calculate LAI for standard output
-    param <- utils::read.table(file.path(gsub(pattern = "/out/",
-                                              replacement = "/run/", x = outdir),
-                                         "sipnet.param"), stringsAsFactors = FALSE)
-    id <- which(param[, 1] == "leafCSpWt")
-    leafC <- 0.48
-    SLA <- 1000 / param[id, 2] #SLA, m2/kgC
-    output[["LAI"]] <- output[["leaf_carbon_content"]] * SLA # LAI
-    output[["fine_root_carbon_content"]] <- sub.sipnet.output$fineRootC   * 0.001  ## fine_root_carbon_content kgC/m2
-    output[["coarse_root_carbon_content"]] <- sub.sipnet.output$coarseRootC * 0.001  ## coarse_root_carbon_content kgC/m2
-    output[["GWBI"]] <- (sub.sipnet.output$woodCreation * 0.001) / 86400 ## kgC/m2/s - this is daily in SIPNET
-    output[["AGB"]] <- (sub.sipnet.output$plantWoodC + sub.sipnet.output$plantLeafC) * 0.001 # Total aboveground biomass kgC/m2
-    output[["time_bounds"]] <- c(rbind(bounds[,1], bounds[,2]))
-    
-    # ******************** Declare netCDF variables ********************#
-    t <- ncdf4::ncdim_def(name = "time",
-                          longname = "time",
-                          units = paste0("days since ", y, "-01-01 00:00:00"),
-                          vals = sub_dates_cf,
-                          calendar = "standard",
-                          unlim = TRUE)
-    lat <- ncdf4::ncdim_def("lat", "degrees_north", vals = as.numeric(sitelat), 
-                            longname = "station_latitude")
-    lon <- ncdf4::ncdim_def("lon", "degrees_east", vals = as.numeric(sitelon), 
-                            longname = "station_longitude")
-    dims <- list(lon = lon, lat = lat, time = t)
-    time_interval <- ncdf4::ncdim_def(name = "hist_interval", 
-                                      longname="history time interval endpoint dimensions",
-                                      vals = 1:2, units="")
-    
-    ## ***** Need to dynamically update the UTC offset here *****
-    
+
+    bounds <- array(data = NA_real_, dim = c(length(sub_dates_cf), 2))
+    bounds[, 1] <- sub_dates_cf
+    bounds[, 2] <- bounds[, 1] + dayfrac
+    bounds <- round(bounds, 4)
+
+    output <- lapply(output_all, function(x) x[row_ids])
+    output[["time_bounds"]] <- c(rbind(bounds[, 1], bounds[, 2]))
+
     for (i in seq_along(output)) {
-      if (length(output[[i]]) == 0)
-        output[[i]] <- rep(-999, length(t$vals))
+      if (length(output[[i]]) == 0) {
+        output[[i]] <- rep(-999, length(sub_dates_cf))
+      }
+      output[[i]][is.na(output[[i]])] <- -999
     }
-    
-    # ******************** Declare netCDF variables ********************#
-    mstmipvar <- PEcAn.utils::mstmipvar
+
+    t <- ncdf4::ncdim_def(
+      name = "time",
+      longname = "time",
+      units = paste0("days since ", y, "-01-01 00:00:00"),
+      vals = sub_dates_cf,
+      calendar = "standard",
+      unlim = TRUE
+    )
+    lat <- ncdf4::ncdim_def("lat", "degrees_north", vals = as.numeric(sitelat), longname = "station_latitude")
+    lon <- ncdf4::ncdim_def("lon", "degrees_east", vals = as.numeric(sitelon), longname = "station_longitude")
+    dims <- list(lon = lon, lat = lat, time = t)
+    time_interval <- ncdf4::ncdim_def(
+      name = "hist_interval",
+      longname = "history time interval endpoint dimensions",
+      vals = 1:2,
+      units = ""
+    )
+
     nc_var <- list(
       "GPP" = PEcAn.utils::to_ncvar("GPP", dims),
       "NPP" = PEcAn.utils::to_ncvar("NPP", dims),
@@ -212,7 +344,7 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
       "AutoResp" = PEcAn.utils::to_ncvar("AutoResp", dims),
       "HeteroResp" = PEcAn.utils::to_ncvar("HeteroResp", dims),
       "SoilResp" = ncdf4::ncvar_def("SoilResp", units = "kg C m-2 s-1", dim = list(lon, lat, t), missval = -999,
-                                    longname = "Soil Respiration"), #need to figure out standard variable for this output
+                                    longname = "Soil Respiration"),
       "NEE" = PEcAn.utils::to_ncvar("NEE", dims),
       "AbvGrndWood" = PEcAn.utils::to_ncvar("AbvGrndWood", dims),
       "leaf_carbon_content" = PEcAn.utils::to_ncvar("leaf_carbon_content", dims),
@@ -232,49 +364,109 @@ model2netcdf.SIPNET <- function(outdir, sitelat, sitelon, start_date, end_date,
                                 longname = "Gross Woody Biomass Increment"),
       "AGB" = ncdf4::ncvar_def("AGB", units = "kg C m-2", dim = list(lon, lat, t), missval = -999,
                                longname = "Total aboveground biomass"),
-      "time_bounds" = ncdf4::ncvar_def(name="time_bounds", units='',
-                                       longname = "history time interval endpoints", dim=list(time_interval,time = t), 
-                                       prec = "double")              
+      "time_bounds" = ncdf4::ncvar_def(name = "time_bounds", units = "",
+                                       longname = "history time interval endpoints", dim = list(time_interval, time = t),
+                                       prec = "double")
     )
-    
-    # ******************** Create netCDF and output variables ********************#
-    ### Output netCDF data
-    if(conflicted & conflict){
-      nc      <- ncdf4::nc_create(file.path(outdir, paste("current", "nc", sep = ".")), nc_var)
-      ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec=NA)
-      for (key in names(nc_var)) {
-        ncdf4::ncvar_put(nc, nc_var[[key]], output[[key]])
+
+    if (file.exists(target_file)) {
+      unlink(target_file)
+    }
+
+    nc <- ncdf4::nc_create(target_file, nc_var)
+    ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec = NA)
+    for (key in names(nc_var)) {
+      ncdf4::ncvar_put(nc, nc_var[[key]], output[[key]])
+    }
+    ncdf4::nc_close(nc)
+
+    invisible(target_file)
+  }
+
+  t_year_loop <- tic()
+  if (!conflict) {
+    years_to_write <- Filter(function(y) {
+      row_ids <- year_index[[as.character(y)]]
+      if (is.null(row_ids) || length(row_ids) == 0) {
+        return(FALSE)
+      }
+      target_file <- file.path(outdir, paste(y, "nc", sep = "."))
+      if (file.exists(target_file) && !overwrite) {
+        return(FALSE)
+      }
+      TRUE
+    }, year_seq)
+
+    if (worker_setting > 1L && length(years_to_write) > 1L) {
+      parallel::mclapply(
+        years_to_write,
+        function(y) {
+          row_ids <- year_index[[as.character(y)]]
+          target_file <- file.path(outdir, paste(y, "nc", sep = "."))
+          write_year_file(y, row_ids, target_file)
+          NULL
+        },
+        mc.cores = min(worker_setting, length(years_to_write))
+      )
+    } else {
+      for (y in years_to_write) {
+        row_ids <- year_index[[as.character(y)]]
+        target_file <- file.path(outdir, paste(y, "nc", sep = "."))
+        write_year_file(y, row_ids, target_file)
+      }
+    }
+  } else {
+    for (y in year_seq) {
+      row_ids <- year_index[[as.character(y)]]
+      if (is.null(row_ids) || length(row_ids) == 0) {
+        next
       }
-      ncdf4::nc_close(nc)
-      
-      #merge nc files of the same year together to enable the assimilation of sub-annual data
-      if(file.exists(file.path(outdir, "previous.nc"))){
-        files <- c(file.path(outdir, "previous.nc"), file.path(outdir, "current.nc"))
-      }else{
-        files <- file.path(outdir, "current.nc")
+
+      target_file <- file.path(outdir, paste(y, "nc", sep = "."))
+      existing_file <- file.exists(target_file)
+
+      if (existing_file && overwrite) {
+        unlink(target_file)
+        existing_file <- FALSE
       }
-      mergeNC(files = files, outfile = file.path(outdir, paste(y, "nc", sep = ".")))
-      #The command "cdo" in mergeNC will automatically rename "time_bounds" to "time_bnds". However, "time_bounds" is used 
-      #in read_restart codes later. So we need to read the new NetCDF file and convert the variable name back. 
-      nc<- ncdf4::nc_open(file.path(outdir, paste(y, "nc", sep = ".")),write=TRUE)
-      nc<-ncdf4::ncvar_rename(nc,"time_bnds","time_bounds")
-      ncdf4::ncatt_put(nc, "time", "bounds","time_bounds", prec=NA)
-      ncdf4::nc_close(nc)
-      unlink(files, recursive = T)
-    }else{
-      nc      <- ncdf4::nc_create(file.path(outdir, paste(y, "nc", sep = ".")), nc_var)
-      ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec=NA)
-      for (i in seq_along(nc_var)) {
-        ncdf4::ncvar_put(nc, nc_var[[i]], output[[i]])
+
+      if (existing_file) {
+        current_file <- tempfile(tmpdir = outdir, pattern = paste0("current_", y, "_"), fileext = ".nc")
+        merged_file <- tempfile(tmpdir = outdir, pattern = paste0("merged_", y, "_"), fileext = ".nc")
+
+        write_year_file(y, row_ids, current_file)
+        mergeNC(files = c(target_file, current_file), outfile = merged_file)
+
+        if (file.exists(merged_file)) {
+          nc <- ncdf4::nc_open(merged_file, write = TRUE)
+          if ("time_bnds" %in% names(nc$var)) {
+            nc <- ncdf4::ncvar_rename(nc, "time_bnds", "time_bounds")
+          }
+          ncdf4::ncatt_put(nc, "time", "bounds", "time_bounds", prec = NA)
+          ncdf4::nc_close(nc)
+
+          if (!file.rename(merged_file, target_file)) {
+            copied <- file.copy(merged_file, target_file, overwrite = TRUE)
+            unlink(merged_file)
+            if (!copied) {
+              PEcAn.logger::logger.error("Failed to replace target file after merge", target_file)
+            }
+          }
+        }
+
+        unlink(current_file)
+      } else {
+        write_year_file(y, row_ids, target_file)
       }
-      ncdf4::nc_close(nc)
     }
-  }  ### End of year loop
-  
-  ## Delete raw output, if requested
+  }
+  toc("yearly_netcdf_conversion", t_year_loop)
+
   if (delete.raw) {
     file.remove(sipnet_out_file)
   }
-} # model2netcdf.SIPNET
+
+  toc("model2netcdf_total", t_total)
+}
 #--------------------------------------------------------------------------------------------------#
 ### EOF
diff --git a/models/sipnet/R/write.configs.SIPNET.R b/models/sipnet/R/write.configs.SIPNET.R
index 3fae8216f18..2e58ae644c4 100755
--- a/models/sipnet/R/write.configs.SIPNET.R
+++ b/models/sipnet/R/write.configs.SIPNET.R
@@ -12,11 +12,55 @@
 ##' @export
 ##' @importFrom rlang %||%
 ##' @author Michael Dietze
+.sipnet_config_cache <- new.env(parent = emptyenv())
+
+.sipnet_profile_enabled <- function() {
+  tolower(trimws(Sys.getenv("PECAN_SIPNET_PROFILE", ""))) %in% c("1", "true", "yes", "on")
+}
+
+.sipnet_verbose_enabled <- function() {
+  tolower(trimws(Sys.getenv("PECAN_SIPNET_VERBOSE", ""))) %in% c("1", "true", "yes", "on")
+}
+
+.sipnet_cache_read_lines <- function(path) {
+  key <- paste0("lines::", normalizePath(path, mustWork = FALSE))
+  if (!exists(key, envir = .sipnet_config_cache, inherits = FALSE)) {
+    assign(key, readLines(con = path, n = -1), envir = .sipnet_config_cache)
+  }
+  get(key, envir = .sipnet_config_cache, inherits = FALSE)
+}
+
+.sipnet_cache_read_table <- function(path) {
+  key <- paste0("table::", normalizePath(path, mustWork = FALSE))
+  if (!exists(key, envir = .sipnet_config_cache, inherits = FALSE)) {
+    assign(key, utils::read.table(path), envir = .sipnet_config_cache)
+  }
+  get(key, envir = .sipnet_config_cache, inherits = FALSE)
+}
+
 write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs = NULL, IC = NULL,
                                 restart = NULL, spinup = NULL) {
+  profile_enabled <- .sipnet_profile_enabled()
+  timings <- list()
+  tic <- function() proc.time()[["elapsed"]]
+  toc <- function(stage, started_at) {
+    if (!profile_enabled) {
+      return(invisible(NULL))
+    }
+    timings[[length(timings) + 1]] <<- data.frame(
+      function_name = "write.config.SIPNET",
+      run_id = as.character(run.id),
+      stage = stage,
+      elapsed_seconds = unname(proc.time()[["elapsed"]] - started_at),
+      stringsAsFactors = FALSE
+    )
+    invisible(NULL)
+  }
+
   ### WRITE sipnet.in
+  t_main <- tic()
   template.in <- system.file("sipnet.in", package = "PEcAn.SIPNET")
-  config.text <- readLines(con = template.in, n = -1)
+  config.text <- .sipnet_cache_read_lines(template.in)
   writeLines(config.text, con = file.path(settings$rundir, run.id, "sipnet.in"))
   
   ### WRITE *.clim
@@ -27,7 +71,9 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       template.clim <- inputs$met$path
     }
   }
-  PEcAn.logger::logger.info(paste0("Writing SIPNET configs with input ", template.clim))
+  if (.sipnet_verbose_enabled()) {
+    PEcAn.logger::logger.info(paste0("Writing SIPNET configs with input ", template.clim))
+  }
   
   # find out where to write run/ouput
   rundir <- file.path(settings$host$rundir, as.character(run.id))
@@ -36,12 +82,50 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     rundir <- file.path(settings$rundir, as.character(run.id))
     outdir <- file.path(settings$modeloutdir, as.character(run.id))
   }
+
+  write_profile_csv <- trimws(Sys.getenv("PECAN_SIPNET_PROFILE_CSV", ""))
+  write_profile_csv <- if (profile_enabled) {
+    if (tolower(write_profile_csv) %in% c("1", "true", "yes", "on")) {
+      file.path(outdir, "sipnet_write.config.profile.csv")
+    } else if (nzchar(write_profile_csv)) {
+      write_profile_csv
+    } else {
+      ""
+    }
+  } else {
+    ""
+  }
+
+  flush_timings <- function() {
+    if (!profile_enabled || length(timings) == 0) {
+      return(invisible(NULL))
+    }
+    timing_df <- do.call(rbind, timings)
+    summary_txt <- apply(timing_df, 1, function(x) {
+      paste0(x[["function_name"]], " [", x[["run_id"]], "] ", x[["stage"]], ": ", signif(as.numeric(x[["elapsed_seconds"]]), 4), " s")
+    })
+    PEcAn.logger::logger.info("SIPNET profiling summary:\n", paste(summary_txt, collapse = "\n"), wrap = FALSE)
+    if (nzchar(write_profile_csv)) {
+      utils::write.table(
+        timing_df,
+        file = write_profile_csv,
+        sep = ",",
+        row.names = FALSE,
+        col.names = !file.exists(write_profile_csv),
+        quote = TRUE,
+        append = file.exists(write_profile_csv)
+      )
+    }
+    invisible(NULL)
+  }
+  on.exit(flush_timings(), add = TRUE)
   
   # create launch script (which will create symlink)
+  t_jobsh <- tic()
   if (!is.null(settings$model$jobtemplate) && file.exists(settings$model$jobtemplate)) {
-    jobsh <- readLines(con = settings$model$jobtemplate, n = -1)
+    jobsh <- .sipnet_cache_read_lines(settings$model$jobtemplate)
   } else {
-    jobsh <- readLines(con = system.file("template.job", package = "PEcAn.SIPNET"), n = -1)
+    jobsh <- .sipnet_cache_read_lines(system.file("template.job", package = "PEcAn.SIPNET"))
   }
   
   # create host specific setttings
@@ -130,6 +214,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   
   writeLines(jobsh, con = file.path(settings$rundir, run.id, "job.sh"))
   Sys.chmod(file.path(settings$rundir, run.id, "job.sh"))
+  toc("write_job_sh", t_jobsh)
   
 
   ### Copy event file
@@ -147,12 +232,22 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
   file.copy(template.paramSpatial, file.path(settings$rundir, run.id, "sipnet.param-spatial"))
   
   ### WRITE *.param
+  t_param <- tic()
   template.param <- system.file("template.param", package = "PEcAn.SIPNET")
   if ("default.param" %in% names(settings$model)) {
     template.param <- settings$model$default.param
   }
   
-  param <- utils::read.table(template.param)
+  param <- .sipnet_cache_read_table(template.param)
+  param_index <- as.list(seq_len(nrow(param)))
+  names(param_index) <- param[, 1]
+  idx <- function(name) {
+    out <- param_index[[name]]
+    if (is.null(out)) {
+      PEcAn.logger::logger.error("Missing parameter", sQuote(name), "in SIPNET template")
+    }
+    out
+  }
   
   #### write run-specific PFT parameters here
   #
@@ -175,14 +270,37 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       "write.config.SIPNET will use the value it sees last."
     )
   }
+  constant.traits <- unlist(defaults[[1]]$constants)
+  constant.names <- names(constant.traits)
+
+  leafphdata <- NULL
+  obs_year_start <- NA_integer_
+  obs_year_end <- NA_integer_
+  if (!is.null(settings$run$inputs$leaf_phenology)) {
+    obs_year_start <- lubridate::year(settings$run$start.date)
+    obs_year_end <- lubridate::year(settings$run$end.date)
+    if (obs_year_start != obs_year_end) {
+      PEcAn.logger::logger.info(
+        "Start.date and end.date are not in the same year.",
+        "Using phenological data from start year only."
+      )
+    }
+    leaf_pheno_path <- settings$run$inputs$leaf_phenology$path
+    if (!is.null(leaf_pheno_path)) {
+      leafphdata <- utils::read.csv(leaf_pheno_path)
+    } else {
+      PEcAn.logger::logger.info("No phenology data were found.",
+        "Please consider running `PEcAn.data.remote::extract_phenology_MODIS`",
+        "to get the parameter file."
+      )
+    }
+  }
+
   for (pft in seq_along(trait.values)) {
     pft.traits <- unlist(trait.values[[pft]])
     pft.trait.names <- names(pft.traits)
     
     ## Append/replace params specified as constants
-    constant.traits <- unlist(defaults[[1]]$constants)
-    constant.names <- names(constant.traits)
-    
     # Replace matches
     for (i in seq_along(constant.traits)) {
       ind <- match(constant.names[i], pft.trait.names)
@@ -206,7 +324,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     if ("leafC" %in% pft.trait.names) {
       leafC <- pft.traits[pft.trait.names == "leafC"] |>
         PEcAn.utils::ud_convert("percent", "1") # percentage to fraction
-      id <- which(param[, 1] == "cFracLeaf")
+      id <- idx("cFracLeaf")
       param[id, 2] <- leafC
     } else {
       leafC <- 0.48 # Fixed value if not available, because it is used in downstream calculations
@@ -214,7 +332,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
 
     # Specific leaf area converted to SLW
     # leafCSpWt [gC/m2 leaf], SLA [m2 leaf/kg leaf], leafC [g C / g leaf]
-    id <- which(param[, 1] == "leafCSpWt")
+    id <- idx("leafCSpWt")
     if ("SLA" %in% pft.trait.names) {
       SLA <- pft.traits[which(pft.trait.names == "SLA")]
       param[id, 2] <- PEcAn.utils::ud_convert(leafC / SLA, "kg/m2", "g/m2")
@@ -225,7 +343,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     # Maximum photosynthesis
     # SIPNET: aMax [nmol CO2 / g   leaf / sec]
     # PEcAn:  Amax [umol CO2 / m^2 leaf / sec]
-    id <- which(param[, 1] == "aMax")
+    id <- idx("aMax")
     SLA_g <- PEcAn.utils::ud_convert(SLA, "1/kg", "1/g") 
     if ("Amax" %in% pft.trait.names) {
       Amax_area <- pft.traits[which(pft.trait.names == "Amax")] # [µmol/m2/s]
@@ -237,34 +355,34 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     
     # Daily fraction of maximum photosynthesis
     if ("AmaxFrac" %in% pft.trait.names) {
-      param[which(param[, 1] == "aMaxFrac"), 2] <- pft.traits[which(pft.trait.names == "AmaxFrac")]
+      param[idx("aMaxFrac"), 2] <- pft.traits[which(pft.trait.names == "AmaxFrac")]
     }
     
     ### Canopy extinction coefficiet (k)
     if ("extinction_coefficient" %in% pft.trait.names) {
-      param[which(param[, 1] == "attenuation"), 2] <- pft.traits[which(pft.trait.names == "extinction_coefficient")]
+      param[idx("attenuation"), 2] <- pft.traits[which(pft.trait.names == "extinction_coefficient")]
     }
     
     # Leaf respiration rate converted to baseFolRespFrac
     if ("leaf_respiration_rate_m2" %in% pft.trait.names) {
       Rd <- pft.traits[which(pft.trait.names == "leaf_respiration_rate_m2")]
-      id <- which(param[, 1] == "baseFolRespFrac")
+      id <- idx("baseFolRespFrac")
       param[id, 2] <- max(min(Rd / Amax_area, 1), 0)
     }
     
     # Low temp threshold for photosynethsis
     if ("Vm_low_temp" %in% pft.trait.names) {
-      param[which(param[, 1] == "psnTMin"), 2] <- pft.traits[which(pft.trait.names == "Vm_low_temp")]
+      param[idx("psnTMin"), 2] <- pft.traits[which(pft.trait.names == "Vm_low_temp")]
     }
     
     # Opt. temp for photosynthesis
     if ("psnTOpt" %in% pft.trait.names) {
-      param[which(param[, 1] == "psnTOpt"), 2] <- pft.traits[which(pft.trait.names == "psnTOpt")]
+      param[idx("psnTOpt"), 2] <- pft.traits[which(pft.trait.names == "psnTOpt")]
     }
     
     # Growth respiration factor (fraction of GPP)
     if ("growth_resp_factor" %in% pft.trait.names) {
-      param[which(param[, 1] == "growthRespFrac"), 2] <- pft.traits[which(pft.trait.names == "growth_resp_factor")]
+      param[idx("growthRespFrac"), 2] <- pft.traits[which(pft.trait.names == "growth_resp_factor")]
     }
     ### !!! NOT YET USED
     #Jmax = NA
@@ -280,7 +398,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     
     # Half saturation of PAR.  PAR at which photosynthesis occurs at 1/2 theoretical maximum (Einsteins * m^-2 ground area * day^-1).
     #if(!is.na(Jmax) & !is.na(alpha)){
-    # param[which(param[,1] == "halfSatPar"),2] = Jmax/(2*alpha)
+    # param[idx("halfSatPar"),2] = Jmax/(2*alpha)
     ### WARNING: this is a very coarse linear approximation and needs improvement *****
     ### Yes, we also need to work on doing a paired query where we have both data together.
     ### Once halfSatPar is calculated, need to remove Jmax and quantum_efficiency from param list so they are not included in SA
@@ -290,45 +408,45 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     # Half saturation of PAR.  PAR at which photosynthesis occurs at 1/2 theoretical maximum (Einsteins * m^-2 ground area * day^-1).
     # Temporary implementation until above is working.
     if ("half_saturation_PAR" %in% pft.trait.names) {
-      param[which(param[, 1] == "halfSatPar"), 2] <- pft.traits[which(pft.trait.names == "half_saturation_PAR")]
+      param[idx("halfSatPar"), 2] <- pft.traits[which(pft.trait.names == "half_saturation_PAR")]
     }
     
     # Ball-berry slomatal slope parameter m
     if ("stomatal_slope.BB" %in% pft.trait.names) {
-      id <- which(param[, 1] == "m_ballBerry")
+      id <- idx("m_ballBerry")
       param[id, 2] <- pft.traits[which(pft.trait.names == "stomatal_slope.BB")]
     }
     
     # Slope of VPD–photosynthesis relationship. dVpd = 1 - dVpdSlope * vpd^dVpdExp
     if ("dVPDSlope" %in% pft.trait.names) {
-      param[which(param[, 1] == "dVpdSlope"), 2] <- pft.traits[which(pft.trait.names == "dVPDSlope")]
+      param[idx("dVpdSlope"), 2] <- pft.traits[which(pft.trait.names == "dVPDSlope")]
     }
     
     # VPD–water use efficiency relationship.  dVpd = 1 - dVpdSlope * vpd^dVpdExp
     if ("dVpdExp" %in% pft.trait.names) {
-      param[which(param[, 1] == "dVpdExp"), 2] <- pft.traits[which(pft.trait.names == "dVpdExp")]
+      param[idx("dVpdExp"), 2] <- pft.traits[which(pft.trait.names == "dVpdExp")]
     }
     
     # Leaf turnover rate average turnover rate of leaves, in fraction per day NOTE: read in as
     # per-year rate!
     if ("leaf_turnover_rate" %in% pft.trait.names) {
-      param[which(param[, 1] == "leafTurnoverRate"), 2] <- pft.traits[which(pft.trait.names == "leaf_turnover_rate")]
+      param[idx("leafTurnoverRate"), 2] <- pft.traits[which(pft.trait.names == "leaf_turnover_rate")]
     }
     
     if ("wueConst" %in% pft.trait.names) {
-      param[which(param[, 1] == "wueConst"), 2] <- pft.traits[which(pft.trait.names == "wueConst")]
+      param[idx("wueConst"), 2] <- pft.traits[which(pft.trait.names == "wueConst")]
     }
     # vegetation respiration Q10.
     if ("veg_respiration_Q10" %in% pft.trait.names) {
-      param[which(param[, 1] == "vegRespQ10"), 2] <- pft.traits[which(pft.trait.names == "veg_respiration_Q10")]
+      param[idx("vegRespQ10"), 2] <- pft.traits[which(pft.trait.names == "veg_respiration_Q10")]
     }
     
     # Base vegetation respiration. vegetation maintenance respiration at 0 degrees C (g C respired * g^-1 plant C * day^-1)
     # NOTE: only counts plant wood C - leaves handled elsewhere (both above and below-ground: assumed for now to have same resp. rate)
     # NOTE: read in as per-year rate!
     if ("stem_respiration_rate" %in% pft.trait.names) {
-      vegRespQ10 <- param[which(param[, 1] == "vegRespQ10"), 2]
-      id <- which(param[, 1] == "baseVegResp")
+      vegRespQ10 <- param[idx("vegRespQ10"), 2]
+      id <- idx("baseVegResp")
       ## Convert from umols CO2 kg s-1 to gC g day-1
       stem_resp_g <- (((pft.traits[which(pft.trait.names == "stem_respiration_rate")]) *
                          (44.0096 / 1e+06) * (12.01 / 44.0096)) / 1000) * 86400
@@ -339,19 +457,19 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     
     # turnover of fine roots (per year rate)
     if ("root_turnover_rate" %in% pft.trait.names) {
-      id <- which(param[, 1] == "fineRootTurnoverRate")
+      id <- idx("fineRootTurnoverRate")
       param[id, 2] <- pft.traits[which(pft.trait.names == "root_turnover_rate")]
     }
     
     # fine root respiration Q10
     if ("fine_root_respiration_Q10" %in% pft.trait.names) {
-      param[which(param[, 1] == "fineRootQ10"), 2] <- pft.traits[which(pft.trait.names == "fine_root_respiration_Q10")]
+      param[idx("fineRootQ10"), 2] <- pft.traits[which(pft.trait.names == "fine_root_respiration_Q10")]
     }
     
     # base respiration rate of fine roots (per year rate)
     if ("root_respiration_rate" %in% pft.trait.names) {
-      fineRootQ10 <- param[which(param[, 1] == "fineRootQ10"), 2]
-      id <- which(param[, 1] == "baseFineRootResp")
+      fineRootQ10 <- param[idx("fineRootQ10"), 2]
+      id <- idx("baseFineRootResp")
       ## Convert from umols CO2 kg s-1 to gC g day-1
       root_resp_rate_g <- (((pft.traits[which(pft.trait.names == "root_respiration_rate")]) *
                               (44.0096/1e+06) * (12.01 / 44.0096)) / 1000) * 86400
@@ -362,7 +480,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     
     # coarse root respiration Q10
     if ("coarse_root_respiration_Q10" %in% pft.trait.names) {
-      param[which(param[, 1] == "coarseRootQ10"), 2] <- pft.traits[which(pft.trait.names == "coarse_root_respiration_Q10")]
+      param[idx("coarseRootQ10"), 2] <- pft.traits[which(pft.trait.names == "coarse_root_respiration_Q10")]
     }
     # WARNING: fineRootAllocation + woodAllocation + leafAllocation isn't supposed to exceed 1
     # see sipnet.c code L2005 :
@@ -383,111 +501,99 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     
     # fineRootAllocation
     if ("root_allocation_fraction" %in% pft.trait.names) {
-      param[which(param[, 1] == "fineRootAllocation"), 2] <- pft.traits[which(pft.trait.names == "root_allocation_fraction")]
+      param[idx("fineRootAllocation"), 2] <- pft.traits[which(pft.trait.names == "root_allocation_fraction")]
     }
     
     # woodAllocation
     if ("wood_allocation_fraction" %in% pft.trait.names) {
-      param[which(param[, 1] == "woodAllocation"), 2] <- pft.traits[which(pft.trait.names == "wood_allocation_fraction")]
+      param[idx("woodAllocation"), 2] <- pft.traits[which(pft.trait.names == "wood_allocation_fraction")]
     }
     
     # leafAllocation
     if ("leaf_allocation_fraction" %in% pft.trait.names) {
-      param[which(param[, 1] == "leafAllocation"), 2] <- pft.traits[which(pft.trait.names == "leaf_allocation_fraction")]
+      param[idx("leafAllocation"), 2] <- pft.traits[which(pft.trait.names == "leaf_allocation_fraction")]
     }
     
     # wood_turnover_rate
     if ("wood_turnover_rate" %in% pft.trait.names) {
-      param[which(param[, 1] == "woodTurnoverRate"), 2] <- pft.traits[which(pft.trait.names == "wood_turnover_rate")]
+      param[idx("woodTurnoverRate"), 2] <- pft.traits[which(pft.trait.names == "wood_turnover_rate")]
     }
     
     ### ----- Soil parameters soil respiration Q10.
     if ("soil_respiration_Q10" %in% pft.trait.names) {
-      param[which(param[, 1] == "soilRespQ10"), 2] <- pft.traits[which(pft.trait.names == "soil_respiration_Q10")]
+      param[idx("soilRespQ10"), 2] <- pft.traits[which(pft.trait.names == "soil_respiration_Q10")]
     }
     # soil respiration rate -- units = 1/year, reference = 0C
     if ("som_respiration_rate" %in% pft.trait.names) {
-      param[which(param[, 1] == "baseSoilResp"), 2] <- pft.traits[which(pft.trait.names == "som_respiration_rate")]
+      param[idx("baseSoilResp"), 2] <- pft.traits[which(pft.trait.names == "som_respiration_rate")]
     }
     
     # litterBreakdownRate
     if ("turn_over_time" %in% pft.trait.names) {
-      id <- which(param[, 1] == "litterBreakdownRate")
+      id <- idx("litterBreakdownRate")
       param[id, 2] <- pft.traits[which(pft.trait.names == "turn_over_time")]
     }
     # frozenSoilEff
     if ("frozenSoilEff" %in% pft.trait.names) {
-      param[which(param[, 1] == "frozenSoilEff"), 2] <- pft.traits[which(pft.trait.names == "frozenSoilEff")]
+      param[idx("frozenSoilEff"), 2] <- pft.traits[which(pft.trait.names == "frozenSoilEff")]
     }
     
     # frozenSoilFolREff
     if ("frozenSoilFolREff" %in% pft.trait.names) {
-      param[which(param[, 1] == "frozenSoilFolREff"), 2] <- pft.traits[which(pft.trait.names == "frozenSoilFolREff")]
+      param[idx("frozenSoilFolREff"), 2] <- pft.traits[which(pft.trait.names == "frozenSoilFolREff")]
     }
     
     # soilWHC
     if ("soilWHC" %in% pft.trait.names) {
-      param[which(param[, 1] == "soilWHC"), 2] <- pft.traits[which(pft.trait.names == "soilWHC")]
+      param[idx("soilWHC"), 2] <- pft.traits[which(pft.trait.names == "soilWHC")]
     }
     # 10/31/2017 IF: these were the two assumptions used in the emulator paper in order to reduce dimensionality
     # These results in improved winter soil respiration values
     # they don't affect anything when the seasonal soil respiration functionality in SIPNET is turned-off
     if(TRUE){
       # assume soil resp Q10 cold == soil resp Q10
-      param[which(param[, 1] == "soilRespQ10Cold"), 2] <- param[which(param[, 1] == "soilRespQ10"), 2]
+      param[idx("soilRespQ10Cold"), 2] <- param[idx("soilRespQ10"), 2]
       # default SIPNET prior of baseSoilRespCold was 1/4th of baseSoilResp
       # assuming they will scale accordingly
-      param[which(param[, 1] == "baseSoilRespCold"), 2] <- param[which(param[, 1] == "baseSoilResp"), 2] * 0.25
+      param[idx("baseSoilRespCold"), 2] <- param[idx("baseSoilResp"), 2] * 0.25
     }
     
     if ("immedEvapFrac" %in% pft.trait.names) {
-      param[which(param[, 1] == "immedEvapFrac"), 2] <- pft.traits[which(pft.trait.names == "immedEvapFrac")]
+      param[idx("immedEvapFrac"), 2] <- pft.traits[which(pft.trait.names == "immedEvapFrac")]
     }
     
     if ("leafWHC" %in% pft.trait.names) {
-      param[which(param[, 1] == "leafPoolDepth"), 2] <- pft.traits[which(pft.trait.names == "leafWHC")]
+      param[idx("leafPoolDepth"), 2] <- pft.traits[which(pft.trait.names == "leafWHC")]
     }
     
     if ("waterRemoveFrac" %in% pft.trait.names) {
-      param[which(param[, 1] == "waterRemoveFrac"), 2] <- pft.traits[which(pft.trait.names == "waterRemoveFrac")]
+      param[idx("waterRemoveFrac"), 2] <- pft.traits[which(pft.trait.names == "waterRemoveFrac")]
     }
     
     if ("fastFlowFrac" %in% pft.trait.names) {
-      param[which(param[, 1] == "fastFlowFrac"), 2] <- pft.traits[which(pft.trait.names == "fastFlowFrac")]
+      param[idx("fastFlowFrac"), 2] <- pft.traits[which(pft.trait.names == "fastFlowFrac")]
     }
     
     if ("rdConst" %in% pft.trait.names) {
-      param[which(param[, 1] == "rdConst"), 2] <- pft.traits[which(pft.trait.names == "rdConst")]
+      param[idx("rdConst"), 2] <- pft.traits[which(pft.trait.names == "rdConst")]
     }
     ### ----- Phenology parameters GDD leaf on
     if ("GDD" %in% pft.trait.names) {
-      param[which(param[, 1] == "gddLeafOn"), 2] <- pft.traits[which(pft.trait.names == "GDD")]
+      param[idx("gddLeafOn"), 2] <- pft.traits[which(pft.trait.names == "GDD")]
     }
     
     # Fraction of leaf fall per year (should be 1 for decid)
     if ("fracLeafFall" %in% pft.trait.names) {
-      param[which(param[, 1] == "fracLeafFall"), 2] <- pft.traits[which(pft.trait.names == "fracLeafFall")]
+      param[idx("fracLeafFall"), 2] <- pft.traits[which(pft.trait.names == "fracLeafFall")]
     }
     
     # Leaf growth.  Amount of C added to the leaf during the greenup period
     if ("leafGrowth" %in% pft.trait.names) {
-      param[which(param[, 1] == "leafGrowth"), 2] <- pft.traits[which(pft.trait.names == "leafGrowth")]
-    }
-
-    #update LeafOnday and LeafOffDay
-    if (!is.null(settings$run$inputs$leaf_phenology)) {
-      obs_year_start <- lubridate::year(settings$run$start.date)
-      obs_year_end <- lubridate::year(settings$run$end.date)
-      if (obs_year_start != obs_year_end) {
-        PEcAn.logger::logger.info(
-          "Start.date and end.date are not in the same year.",
-          "Using phenological data from start year only."
-        )
-      }
-      leaf_pheno_path <- settings$run$inputs$leaf_phenology$path
-      if (!is.null(leaf_pheno_path)) {
-        ##read data
-        leafphdata <- utils::read.csv(leaf_pheno_path) #leaf phenology data starting from 2001-01-01 to current
+      param[idx("leafGrowth"), 2] <- pft.traits[which(pft.trait.names == "leafGrowth")]
+    }
+
+    # update LeafOnday and LeafOffDay
+    if (!is.null(leafphdata)) {
         leafOnDay <- leafphdata$leafonday[leafphdata$year == obs_year_start
                                           & leafphdata$site_id == settings$run$site$id]
         leafOffDay <- leafphdata$leafoffday[leafphdata$year == obs_year_start
@@ -503,7 +609,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
             PEcAn.logger::logger.info(paste("Missing leafOnDay for current year. Using site mean:", leafOnDay))
           } else {
             # 2. If no site history exists, fall back to parameter file
-            leafOnDay <- param[which(param[, 1] == "leafOnDay"), 2]
+            leafOnDay <- param[idx("leafOnDay"), 2]
             PEcAn.logger::logger.warn("Missing leafOnDay and no site history. Using parameter file default.")
           }
         }
@@ -518,7 +624,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
             PEcAn.logger::logger.info(paste("Missing leafOffDay for current year. Using site mean:", leafOffDay))
           } else {
             # 2. If no site history exists, fall back to parameter file
-            leafOffDay <- param[which(param[, 1] == "leafOffDay"), 2]
+            leafOffDay <- param[idx("leafOffDay"), 2]
             PEcAn.logger::logger.warn("Missing leafOffDay and no site history. Using parameter file default.")
           }
         }
@@ -526,18 +632,15 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
         # when we have Leaf off date larger than leaf on date.
         # Otherwise the phenology will not be used.
         if (leafOffDay > leafOnDay) {
-          param[which(param[, 1] == "leafOnDay"), 2] <- leafOnDay
-          param[which(param[, 1] == "leafOffDay"), 2] <- leafOffDay
+          param[idx("leafOnDay"), 2] <- leafOnDay
+          param[idx("leafOffDay"), 2] <- leafOffDay
         }
-      } else {
-        PEcAn.logger::logger.info("No phenology data were found.",
-          "Please consider running `PEcAn.data.remote::extract_phenology_MODIS`",
-          "to get the parameter file."
-        )
-      }
     }
   } ## end loop over PFTS
+  toc("update_params_from_traits", t_param)
   ####### end parameter update
+
+  t_state_inputs <- tic()
   #working on reading soil file
   if (length(settings$run$inputs$soil_physics$path) > 0) {
     template.soil_physics <- settings$run$inputs$soil_physics$path  ## read from settings
@@ -574,7 +677,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
           soilWHC_total <- sum(unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])*thickness)
           if (thickness[1]<=10) {
             #LitterWHC in cm, assuming the litter depth is within the top 10 cm
-            param[which(param[, 1] == "litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*thickness[1]
+            param[idx("litterWHC"), 2] <- unlist(soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"])[1]*thickness[1]
           }
         } else {
           #if no depth/thickness is provided
@@ -582,11 +685,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
           thickness <- 100 #assume the default soil depth is the plant rooting depth of 100 cm, or use the user-specified value
           soilWHC_total <- soil_IC_list$vals["volume_fraction_of_water_in_soil_at_saturation"]*thickness
         }
-        param[which(param[, 1] == "soilWHC"), 2] <- soilWHC_total
+        param[idx("soilWHC"), 2] <- soilWHC_total
       }
       if ("soil_hydraulic_conductivity_at_saturation" %in% names(soil_IC_list$vals)) {
          #litwaterDrainrate in cm/day
-         param[which(param[, 1] == "litWaterDrainRate"), 2] <- PEcAn.utils::ud_convert(unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1], "m s-1", "cm day-1")
+         param[idx("litWaterDrainRate"), 2] <- PEcAn.utils::ud_convert(unlist(soil_IC_list$vals["soil_hydraulic_conductivity_at_saturation"])[1], "m s-1", "cm day-1")
        }
     }
   }
@@ -618,42 +721,42 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
             )
           )
         }
-      param[which(param[, 1] == "plantWoodInit"),  2] <- wood_total_C
-      param[which(param[, 1] == "coarseRootFrac"), 2] <- IC$coarseRootFrac
-      param[which(param[, 1] == "fineRootFrac"),   2] <- IC$fineRootFrac
+      param[idx("plantWoodInit"),  2] <- wood_total_C
+      param[idx("coarseRootFrac"), 2] <- IC$coarseRootFrac
+      param[idx("fineRootFrac"),   2] <- IC$fineRootFrac
     }
     ## laiInit m2/m2
     if ("lai" %in% ic.names) {
-      param[which(param[, 1] == "laiInit"), 2] <- IC$lai
+      param[idx("laiInit"), 2] <- IC$lai
     }
     ## litterInit gC/m2
     if ("litter_carbon_content" %in% ic.names) {
-      param[which(param[, 1] == "litterInit"), 2] <- IC$litter_carbon_content
+      param[idx("litterInit"), 2] <- IC$litter_carbon_content
     }
     ## soilInit gC/m2
     if ("soil" %in% ic.names) {
-      param[which(param[, 1] == "soilInit"), 2] <- IC$soil
+      param[idx("soilInit"), 2] <- IC$soil
     }
     ## litterWFracInit fraction
     if ("litter_mass_content_of_water" %in% ic.names) {
       #here we use litterWaterContent/litterWHC to calculate the litterWFracInit
-      param[which(param[, 1] == "litterWFracInit"), 2] <- IC$litter_mass_content_of_water/(param[which(param[, 1] == "litterWHC"), 2]*10)
+      param[idx("litterWFracInit"), 2] <- IC$litter_mass_content_of_water/(param[idx("litterWHC"), 2]*10)
     }
     ## soilWater IC$soilWater is in kg/m2, and soilWHC is in cm
     if ("soilWater" %in% ic.names) {
-      param[which(param[, 1] == "soilWFracInit"), 2] <- IC$soilWater/(param[which(param[, 1] == "soilWHC"), 2]*10)
+      param[idx("soilWFracInit"), 2] <- IC$soilWater/(param[idx("soilWHC"), 2]*10)
     }
     ## soilWFracInit fraction
     if ("soilWFrac" %in% ic.names) {
-      param[which(param[, 1] == "soilWFracInit"), 2] <- IC$soilWFrac
+      param[idx("soilWFracInit"), 2] <- IC$soilWFrac
     }
     ## snowInit cm water equivalent
     if ("SWE" %in% ic.names) {
-      param[which(param[, 1] == "snowInit"), 2] <- IC$SWE
+      param[idx("snowInit"), 2] <- IC$SWE
     }
     ## microbeInit mgC/g soil
     if ("microbe" %in% ic.names) {
-      param[which(param[, 1] == "microbeInit"), 2] <- IC$microbe
+      param[idx("microbeInit"), 2] <- IC$microbe
     }
 
   } else if (length(settings$run$inputs$poolinitcond$path) > 0) {
@@ -688,15 +791,15 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
 
       ## plantWoodInit gC/m2
       if ("wood" %in% names(IC.pools)) {
-        fineRootFrac <- param[which(param[,1] == "fineRootFrac"),2]
-        coarseRootFrac <- param[which(param[,1] == "coarseRootFrac"),2]
+        fineRootFrac <- param[idx("fineRootFrac"),2]
+        coarseRootFrac <- param[idx("coarseRootFrac"),2]
         # accounts for the fact that SIPNET take plantWoodInit as all woods (including roots).
-        param[which(param[, 1] == "plantWoodInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$wood, "kg m-2", "g m-2")/(1-fineRootFrac-coarseRootFrac)
+        param[idx("plantWoodInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$wood, "kg m-2", "g m-2")/(1-fineRootFrac-coarseRootFrac)
       }
       ## laiInit m2/m2
       lai <- IC.pools$LAI
       if (!is.na(lai) && is.numeric(lai)) {
-        param[param[, 1] == "laiInit", 2] <- lai
+        param[idx("laiInit"), 2] <- lai
       }
 
       # Sipnet always starts from initial LAI whether day 0 is in or out of the
@@ -708,40 +811,40 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       # - the PFT sets leafOnDay/leafOffday as traits.
       # So unless you set something different, it's probably using DOY 144/285
       # ==> leaves are on from late May through mid-October.
-      is_deciduous_pft <- isTRUE(param[param[, 1] == "fracLeafFall", 2] > 0.5)
+      is_deciduous_pft <- isTRUE(param[idx("fracLeafFall"), 2] > 0.5)
       start_day <- lubridate::yday(settings$run$start.date)
       starts_with_leaves <- (
-        start_day >= param[param[, 1] == "leafOnDay", 2]
-        && start_day <= param[param[, 1] == "leafOffDay", 2]
+        start_day >= param[idx("leafOnDay"), 2]
+        && start_day <= param[idx("leafOffDay"), 2]
       )
       if (is_deciduous_pft && !starts_with_leaves) {
         # Note that this doesn't adjust for winter LAI of evergreens!
         # Could consider using LAI*fracLeafFall,
         # But that strongly assumes that IC LAI is both (1) reported at
         # season peak and not (2) adjusted by any earlier step (i.e. SDA).
-        param[param[, 1] == "laiInit", 2] <- 0
+        param[idx("laiInit"), 2] <- 0
       }
 
       ## neeInit gC/m2
       if (ic_has_ncvars[["nee"]]) {
         nee <- ncdf4::ncvar_get(IC.nc, "nee")
         if (!is.na(nee) && is.numeric(nee)) {
-          param[param[, 1] == "neeInit", 2] <- nee
+          param[idx("neeInit"), 2] <- nee
         }
       }
       ## litterInit gC/m2
       if ("litter" %in% names(IC.pools)) {
-        param[param[, 1] == "litterInit", 2] <- PEcAn.utils::ud_convert(IC.pools$litter, "g m-2", "g m-2") # BETY: kgC m-2
+        param[idx("litterInit"), 2] <- PEcAn.utils::ud_convert(IC.pools$litter, "g m-2", "g m-2") # BETY: kgC m-2
       }
       ## soilInit gC/m2
       if ("soil" %in% names(IC.pools)) {
-        param[param[, 1] == "soilInit", 2] <- PEcAn.utils::ud_convert(sum(IC.pools$soil), "kg m-2", "g m-2") # BETY: kgC m-2
+        param[idx("soilInit"), 2] <- PEcAn.utils::ud_convert(sum(IC.pools$soil), "kg m-2", "g m-2") # BETY: kgC m-2
       }
       ## soilWFracInit fraction
       if (ic_has_ncvars[["SoilMoistFrac"]]) {
         soilWFrac <- ncdf4::ncvar_get(IC.nc, "SoilMoistFrac")
         if (!is.na(soilWFrac) && is.numeric(soilWFrac)) {
-          param[param[, 1] == "soilWFracInit", 2] <- sum(soilWFrac) / 100
+          param[idx("soilWFracInit"), 2] <- sum(soilWFrac) / 100
           ## litterWFracInit fraction
           litterWFrac <- soilWFrac
         }
@@ -752,27 +855,27 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       if (ic_has_ncvars[["SWE"]]) {
         snow <- ncdf4::ncvar_get(IC.nc, "SWE")
         if (!is.na(snow) && is.numeric(snow)) {
-          param[param[, 1] == "snowInit", 2] <- PEcAn.utils::ud_convert(snow, "kg m-2", "g cm-2")  # BETY: kg m-2
+          param[idx("snowInit"), 2] <- PEcAn.utils::ud_convert(snow, "kg m-2", "g cm-2")  # BETY: kg m-2
         }
       }
       ## leafOnDay
       if (ic_has_ncvars[["date_of_budburst"]]) {
         leafOnDay <- ncdf4::ncvar_get(IC.nc, "date_of_budburst")
         if (!is.na(leafOnDay) && is.numeric(leafOnDay)) {
-          param[param[, 1] == "leafOnDay", 2] <- leafOnDay
+          param[idx("leafOnDay"), 2] <- leafOnDay
         }
       }
       ## leafOffDay
       if (ic_has_ncvars[["date_of_senescence"]]) {
         leafOffDay <- ncdf4::ncvar_get(IC.nc, "date_of_senescence")
         if (!is.na(leafOffDay) && is.numeric(leafOffDay)) {
-          param[param[, 1] == "leafOffDay", 2] <- leafOffDay
+          param[idx("leafOffDay"), 2] <- leafOffDay
         }
       }
       if (ic_has_ncvars[["Microbial Biomass C"]]) {
         microbe <- ncdf4::ncvar_get(IC.nc, "Microbial Biomass C")
         if (!is.na(microbe) && is.numeric(microbe)) {
-          param[param[, 1] == "microbeInit", 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") #BETY: mg microbial C kg-1 soil
+          param[idx("microbeInit"), 2] <- PEcAn.utils::ud_convert(microbe, "mg kg-1", "mg g-1") #BETY: mg microbial C kg-1 soil
         }
       }
 
@@ -791,10 +894,11 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
       soil.path <- settings$run$inputs$soilmoisture$path
       soilWFrac <- ncdf4::ncvar_get(ncdf4::nc_open(soil.path), varid = "mass_fraction_of_unfrozen_water_in_soil_moisture")
 
-      param[which(param[, 1] == "soilWFracInit"), 2] <- soilWFrac
+      param[idx("soilWFracInit"), 2] <- soilWFrac
     }
 
   }
+  toc("apply_state_and_soil_inputs", t_state_inputs)
   if (file.exists(file.path(settings$rundir, run.id, "sipnet.param"))) {
     file.rename(
       file.path(settings$rundir, run.id, "sipnet.param"),
@@ -814,6 +918,7 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
     col.names = FALSE,
     quote = FALSE
   )
+  toc("write_config_total", t_main)
 } # write.config.SIPNET
 
 
diff --git a/models/sipnet/README.md b/models/sipnet/README.md
index 8fee3a08770..d4affdc31b3 100644
--- a/models/sipnet/README.md
+++ b/models/sipnet/README.md
@@ -37,3 +37,15 @@ This is a basic example which shows you how to solve a common problem:
 library(PEcAn.SIPNET)
 ## basic example code
 ```
+
+## Runtime / Performance Notes
+
+- `data.table` is a required dependency. `model2netcdf.SIPNET` uses `data.table::fread` for `sipnet.out` parsing.
+- `model2netcdf.SIPNET` supports two `sipnet.out` header styles:
+  - header on first line
+  - `Notes:` preamble line before header
+- Internal NetCDF worker default is `max(1, parallel::detectCores() - 1)`.
+  - Override with `PECAN_SIPNET_NC_WORKERS=<N>`.
+- `conflict=TRUE` requires `cdo` in `PATH`; function fails fast if `cdo` is unavailable.
+- `PECAN_SIPNET_PROFILE=1` enables timing summaries (optional CSV via `PECAN_SIPNET_PROFILE_CSV`).
+- `PECAN_SIPNET_VERBOSE=1` enables per-run config input logging in `write.config.SIPNET`.

From 9cd0d8fd3e313876bf60ed13dfe9afb3502f6b4c Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:01:28 -0700
Subject: [PATCH 02/25] replace stop with pecan logger in run.write.configs

---
 base/workflow/R/run.write.configs.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R
index 31e0a388d29..83d1913f0f4 100644
--- a/base/workflow/R/run.write.configs.R
+++ b/base/workflow/R/run.write.configs.R
@@ -32,9 +32,9 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
 
   # Validate that input_design matches ensemble.size
   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 for each run."
+      ensemble.size, ". The design matrix must have exactly one row for each run."
     )
   }
                               

From 61b413b29bf0f8a50c64406e79df5bf021ea49d8 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:04:06 -0700
Subject: [PATCH 03/25] Improve input design validation and logging in
 runModule.run.write.configs

---
 base/workflow/R/runModule.run.write.configs.R | 40 +++++++++++--------
 1 file changed, 23 insertions(+), 17 deletions(-)

diff --git a/base/workflow/R/runModule.run.write.configs.R b/base/workflow/R/runModule.run.write.configs.R
index ae6ef91a90c..cd8aa9b4b16 100644
--- a/base/workflow/R/runModule.run.write.configs.R
+++ b/base/workflow/R/runModule.run.write.configs.R
@@ -7,9 +7,13 @@
 #'   rows from `trait.samples`/`ensemble.samples` plus optional columns named for
 #'   `settings$run$inputs` tags (e.g. `met`, `soil`) with index (i.e., row number)
 #'   into each input's `path` list. Provide at least one row per planned run
-#'   (median + all SA members and/or `ensemble.size`). Usually generated by
+#'   (median + all SA members and/or `ensemble.size`). When supplied,
+#'   `input_design` defines the runs to execute and the effective run count is
+#'   `nrow(input_design)`. In that case, `settings$ensemble$size` is ignored as
+#'   an input and updated to match the design row count. If `input_design` is
+#'   NULL, `settings$ensemble$size` is used to generate a design via
+#'   `generate_joint_ensemble_design()`. Usually generated by
 #'   `generate_joint_ensemble_design()` but custom designs may be supplied.
-#'  If NULL, `generate_joint_ensemble_design()` will be called internally.
 #' @return A modified settings object, invisibly
 #' @importFrom dplyr %>%
 #' @export
@@ -30,13 +34,14 @@ runModule.run.write.configs <- function(settings,
         ensemble_size = ensemble_size
       )
       input_design <- design_result$X
+    } else if (!is.null(settings$ensemble$size) &&
+               settings$ensemble$size != nrow(input_design)) {
+      PEcAn.logger::logger.warn(
+        "Ignoring settings$ensemble$size = ", settings$ensemble$size,
+        "; using nrow(input_design) = ", nrow(input_design), "."
+      )
     }
-
-    # Validate design matrix size for MultiSettings
-    if (!is.null(settings$ensemble$size) && nrow(input_design) != settings$ensemble$size) {
-      PEcAn.logger::logger.severe("Input_design has", nrow(input_design), "rows but settings$ensemble$size is",
-           settings$ensemble$size, ". Design matrix must have exactly one row per run.")
-    }
+    settings$ensemble$size <- nrow(input_design)
 
     return(PEcAn.settings::papply(settings,
                                   runModule.run.write.configs,
@@ -54,16 +59,15 @@ runModule.run.write.configs <- function(settings,
         ensemble_size = ensemble_size
       )
       input_design <- design_result$X
+    } else if (!is.null(settings$ensemble$size) &&
+               settings$ensemble$size != nrow(input_design)) {
+      PEcAn.logger::logger.warn(
+        "Ignoring settings$ensemble$size = ", settings$ensemble$size,
+        "; using nrow(input_design) = ", nrow(input_design), "."
+      )
     }
-
-    # Validate design matrix size for Settings
-    if (!is.null(settings$ensemble$size) && nrow(input_design) != settings$ensemble$size) {
-      PEcAn.logger::logger.severe("Input_design has", nrow(input_design), "rows but settings$ensemble$size is",
-           settings$ensemble$size, ". Design matrix must have exactly one row per run.")
-    }
-
     ensemble_size <- nrow(input_design)
-
+    settings$ensemble$size <- ensemble_size
 
     # check to see if there are posterior.files tags under pft
     posterior.files <- settings$pfts %>%
@@ -78,6 +82,8 @@ runModule.run.write.configs <- function(settings,
       input_design = input_design
     ))
   } else {
-    stop("runModule.run.write.configs only works with Settings or MultiSettings")
+    PEcAn.logger::logger.error(
+      "runModule.run.write.configs only works with Settings or MultiSettings"
+    )
   }
 }

From 7e6450031f668ab688e1d08eadf3cd4c81708a80 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:11:07 -0700
Subject: [PATCH 04/25] Update input design documentation and change
 sensitivity import to sensobol

---
 base/workflow/R/runModule.run.write.configs.R | 18 +++++++-----------
 modules/uncertainty/DESCRIPTION               |  2 +-
 2 files changed, 8 insertions(+), 12 deletions(-)

diff --git a/base/workflow/R/runModule.run.write.configs.R b/base/workflow/R/runModule.run.write.configs.R
index cd8aa9b4b16..839fe5a7ea9 100644
--- a/base/workflow/R/runModule.run.write.configs.R
+++ b/base/workflow/R/runModule.run.write.configs.R
@@ -3,17 +3,13 @@
 #' @param settings a PEcAn Settings or MultiSettings object
 #' @param overwrite logical: Replace config files if they already exist?
 #' @param input_design data.frame design matrix linking parameter draws and any
-#'   sampled inputs across runs. Include a `param` column whose values select
-#'   rows from `trait.samples`/`ensemble.samples` plus optional columns named for
-#'   `settings$run$inputs` tags (e.g. `met`, `soil`) with index (i.e., row number)
-#'   into each input's `path` list. Provide at least one row per planned run
-#'   (median + all SA members and/or `ensemble.size`). When supplied,
-#'   `input_design` defines the runs to execute and the effective run count is
-#'   `nrow(input_design)`. In that case, `settings$ensemble$size` is ignored as
-#'   an input and updated to match the design row count. If `input_design` is
-#'   NULL, `settings$ensemble$size` is used to generate a design via
-#'   `generate_joint_ensemble_design()`. Usually generated by
-#'   `generate_joint_ensemble_design()` but custom designs may be supplied.
+#'   sampled inputs across runs. Must include a `param` column selecting rows
+#'   from `trait.samples`/`ensemble.samples`; may also include columns named for
+#'   `settings$run$inputs` tags (e.g. `met`, `soil`) giving row indices into
+#'   each input's `path` list. If supplied, `input_design` defines the runs to
+#'   execute and the run count is `nrow(input_design)`, overriding
+#'   `settings$ensemble$size`. If NULL, the design is generated from
+#'   `settings$ensemble$size` using `generate_joint_ensemble_design()`.
 #' @return A modified settings object, invisibly
 #' @importFrom dplyr %>%
 #' @export
diff --git a/modules/uncertainty/DESCRIPTION b/modules/uncertainty/DESCRIPTION
index 6a32d353c6d..3ffda9157e9 100644
--- a/modules/uncertainty/DESCRIPTION
+++ b/modules/uncertainty/DESCRIPTION
@@ -42,7 +42,7 @@ Imports:
     purrr,
     randtoolbox,
     rlang,
-    sensitivity
+    sensobol
 Suggests:
     testthat (>= 1.0.2),
     mockery

From 3da3350adbf75b19d4b21115fc752259f65f074b Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:11:54 -0700
Subject: [PATCH 05/25] test that runModule.run.write.configs uses input_design
 row count

---
 .../testthat/test.run.write.configs.sobol.R   | 29 +++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 base/workflow/tests/testthat/test.run.write.configs.sobol.R

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..264b135da21
--- /dev/null
+++ b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
@@ -0,0 +1,29 @@
+library(testthat)
+
+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))
+
+  mockery::stub(
+    runModule.run.write.configs,
+    "PEcAn.workflow::run.write.configs",
+    function(settings, ensemble.size, input_design, write, posterior.files, overwrite) {
+      list(settings = settings, ensemble.size = ensemble.size, input_design = input_design)
+    }
+  )
+
+  result <- runModule.run.write.configs(
+    settings = settings,
+    input_design = input_design
+  )
+
+  expect_equal(result$ensemble.size, 5)
+  expect_equal(result$settings$ensemble$size, 5)
+})

From 1b1b3ac1ff17e3fbbefc1c201949de40c2f5cf29 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:27:36 -0700
Subject: [PATCH 06/25] Switch Sobol postprocessing to sensobol and
 standardized outputs

Replace the old sensitivity-based Sobol postprocessing with sensobol and compute indices from PEcAn standardized ensemble.output.*.Rdata files instead of reconstructing responses from runs.txt and raw per-run outputs.

This keeps Sobol analysis aligned with the normal PEcAn workflow and output conventions while saving reproducible Sobol design and index artifacts.
---
 modules/uncertainty/R/compute_sobol_indices.R | 153 ++++++++++++++----
 1 file changed, 118 insertions(+), 35 deletions(-)

diff --git a/modules/uncertainty/R/compute_sobol_indices.R b/modules/uncertainty/R/compute_sobol_indices.R
index d565a771f63..b6db1864712 100644
--- a/modules/uncertainty/R/compute_sobol_indices.R
+++ b/modules/uncertainty/R/compute_sobol_indices.R
@@ -1,45 +1,128 @@
 #' 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 `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
+#' parameter alone, while total-order indices summarize the full contribution of
+#' that parameter including its interactions with other parameters. PEcAn uses
+#' `sensobol` for these variance-based estimators; readers wanting estimator
+#' details or broader background should refer to Puy et al. (2022).
 #'
-#' @return           sobol_obj 
-#' .
+#' @param outdir PEcAn run output directory containing `ensemble.output.*.Rdata`
+#'   files.
+#' @param sobol_obj object produced by
+#'   `PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)`.
+#' @param var Variable name to summarize (default `"GPP"`).
+#' @param stat_fun Summary statistic applied to `var`. Retained for backwards
+#'   compatibility; standardized ensemble outputs are already scalar summaries.
+#'
+#' @return A tibble of Sobol first-order and total-order indices with attached
+#'   factor metadata.
+#' @references 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))
+  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)
 }

From f4fe4469bfdfb5e87d01a7ab64dec5a58fe68955 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:44:23 -0700
Subject: [PATCH 07/25] add new Sobol design tests in test.sobol.R

---
 .../testthat/test.run.write.configs.sobol.R   |   2 -
 .../uncertainty/tests/testthat/test.sobol.R   | 127 ++++++++++++++++++
 .../tests/testthat/test_ensemble.R            |   1 -
 3 files changed, 127 insertions(+), 3 deletions(-)
 create mode 100644 modules/uncertainty/tests/testthat/test.sobol.R

diff --git a/base/workflow/tests/testthat/test.run.write.configs.sobol.R b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
index 264b135da21..a03b9260820 100644
--- a/base/workflow/tests/testthat/test.run.write.configs.sobol.R
+++ b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
@@ -1,5 +1,3 @@
-library(testthat)
-
 test_that("runModule.run.write.configs uses input_design row count", {
   settings <- PEcAn.settings::Settings(
     ensemble = list(
diff --git a/modules/uncertainty/tests/testthat/test.sobol.R b/modules/uncertainty/tests/testthat/test.sobol.R
new file mode 100644
index 00000000000..e0e6af35ca2
--- /dev/null
+++ b/modules/uncertainty/tests/testthat/test.sobol.R
@@ -0,0 +1,127 @@
+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", parent = "met")
+      )
+    )
+  )
+}
+
+test_that("Sobol design expands to N * (k + 2) rows with metadata", {
+  withr::with_tempdir({
+    settings <- make_sobol_settings(getwd())
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 4,
+      sobol = TRUE
+    )
+
+    expect_equal(nrow(result$X), 20)
+    expect_equal(result$N, 4)
+    expect_identical(result$params, c("param", "met", "poolinitcond"))
+    expect_identical(result$backend, "sensobol")
+    expect_identical(result$matrices, c("A", "B", "AB"))
+    expect_identical(result$first, "saltelli")
+    expect_identical(result$total, "jansen")
+    expect_true(all(c("param", "met", "poolinitcond", "events") %in% names(result$X)))
+    expect_equal(result$X$events, result$X$met)
+    expect_true(all(result$X$param >= 1))
+    expect_true(all(result$X$param <= 20))
+    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"))
+  })
+})
+
+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("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 8197c66d29a..403cc5dc9d9 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

From 55ad31d2579fdc2e29f6859cc29d22a5dcd1cba5 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 15:47:07 -0700
Subject: [PATCH 08/25] Replace implicit Sobol matrix split with explicit
 type-level design generation

---
 .../R/generate_joint_ensemble_design.R        | 169 ++++++++++++++----
 1 file changed, 137 insertions(+), 32 deletions(-)

diff --git a/modules/uncertainty/R/generate_joint_ensemble_design.R b/modules/uncertainty/R/generate_joint_ensemble_design.R
index 02c9289871f..ed7da584791 100644
--- a/modules/uncertainty/R/generate_joint_ensemble_design.R
+++ b/modules/uncertainty/R/generate_joint_ensemble_design.R
@@ -4,24 +4,49 @@
 #' are shared across sites to ensure consistent parameter sampling.
 #'
 #' @param settings A PEcAn settings object containing ensemble configuration
-#' @param ensemble_size Integer specifying the number of ensemble members
-#' Since the `input_design` will only be generated once for the entire model run,
-#' the only situation, where we might want to recycle the existing `ensemble_samples`,
-#' is when we split and submit the larger SDA runs (e.g., 8,000 sites) into 
-#' smaller SDA experiments (e.g., 100 sites per job), where we want to keep using 
-#' the same parameters rather than creating new parameters for each job.
-#' @param sobol for activating sobol
-#' @return  A list containing ensemble samples and indices
-#'   If `sobol = TRUE`, the list will be a `sensitivity::soboljansen()` 
-#'   result and will contain the components documented therein.
+#' @param ensemble_size Integer specifying the number of ensemble members.
+#'   When `sobol = TRUE`, this is the Sobol base sample size `N`, not the
+#'   expanded number of model runs.
+#' @param sobol Logical, generate a variance-based Sobol design using
+#'   `sensobol`.
+#' @return A list with component `X`, a data frame design matrix describing
+#'   PEcAn parameter and sampled-input indices. If `sobol = TRUE`, the list
+#'   also includes the metadata needed by `compute_sobol_indices()`.
 #' @export
 
+.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)
+  }
+
+  first_pft <- samples$trait.samples[[1]]
+  if (is.null(first_pft) || length(first_pft) == 0) {
+    return(0L)
+  }
+
+  first_trait <- first_pft[[1]]
+  if (is.null(first_trait)) {
+    return(0L)
+  }
+
+  return(as.integer(length(first_trait)))
+}
+
+.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 <- 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()
@@ -35,6 +60,99 @@ generate_joint_ensemble_design <- function(settings,
   ]
   samp.ordered <- samp[c(order, names(samp)[!(names(samp) %in% order)])]
 
+  if (sobol) {
+    sobol_factors <- c(
+      "param",
+      names(samp.ordered)[
+        names(samp.ordered) != "parameters" &
+          vapply(
+            samp.ordered,
+            function(x) is.null(x$parent),
+            logical(1)
+          )
+      ]
+    )
+
+    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)
+
+    sobol_indices <- list()
+    sobol_indices[["param"]] <- .map_sobol_to_indices(
+      sobol_design[["param"]],
+      total_runs
+    )
+
+    for (input_tag in setdiff(sobol_factors, "param")) {
+      input_paths <- settings$run$inputs[[tolower(input_tag)]]$path
+      sobol_indices[[input_tag]] <- .map_sobol_to_indices(
+        sobol_design[[input_tag]],
+        length(input_paths)
+      )
+    }
+
+    for (i in seq_along(samp.ordered)) {
+      input_tag <- names(samp.ordered)[i]
+
+      if (identical(input_tag, "parameters")) {
+        next
+      }
+
+      parent_name <- samp.ordered[[i]]$parent
+      if (!is.null(parent_name)) {
+        input_result <- PEcAn.uncertainty::input.ens.gen(
+          settings = settings,
+          ensemble_size = total_runs,
+          input = input_tag,
+          method = samp.ordered[[i]]$method,
+          parent_ids = sampled_inputs[[parent_name]]
+        )
+        sampled_inputs[[input_tag]] <- input_result
+        design_list[[input_tag]] <- input_result$ids
+      } else if (input_tag %in% names(sobol_indices)) {
+        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 <- data.frame(design_list)
+
+    factor_metadata <- data.frame(
+      factor = sobol_factors,
+      source_type = sobol_factors,
+      source_tag = ifelse(sobol_factors == "param", NA_character_, sobol_factors),
+      stringsAsFactors = FALSE
+    )
+
+    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
+    ))
+  }
+
   # loop over inputs.
   for (i in seq_along(samp.ordered)) {
     input_tag <- names(samp.ordered)[i]
@@ -54,33 +172,20 @@ generate_joint_ensemble_design <- function(settings,
       parent_ids = parent_ids
     )
 
-    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.
   return(list(X = design_matrix))
-}
\ No newline at end of file
+}

From f29cdef0ec2a595c0baf8086a41f0c029d824378 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:37:03 -0700
Subject: [PATCH 09/25] Add error handling for trait sample indices in
 run.write.configs - validate input_design$param against trait sample bank

---
 base/workflow/R/run.write.configs.R | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R
index 83d1913f0f4..eedb35eae9b 100644
--- a/base/workflow/R/run.write.configs.R
+++ b/base/workflow/R/run.write.configs.R
@@ -119,6 +119,27 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
     load(samples.file, envir = samples) ## loads ensemble.samples, trait.samples, sa.samples, runs.samples, env.samples
     trait.samples <- samples$trait.samples
     trait_sample_indices <- input_design[["param"]]
+    if (is.null(trait_sample_indices)) {
+      PEcAn.logger::logger.error(
+        "input_design must include a `param` column selecting rows from trait.samples"
+      )
+    }
+    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"
+      )
+    }
+    first_pft <- trait.samples[[1]]
+    first_trait <- first_pft[[1]]
+    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")
+    }
+    if (is.null(first_trait) || any(trait_sample_indices > length(first_trait))) {
+      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]]

From a759c7c4f67dcba76384859874b37f1119491471 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:37:49 -0700
Subject: [PATCH 10/25] Update documentation for compute_sobol_indices: remove
 stat_fun parameter and clarify references

---
 modules/uncertainty/R/compute_sobol_indices.R | 15 +++----
 .../uncertainty/man/compute_sobol_indices.Rd  | 39 +++++++++++++------
 2 files changed, 36 insertions(+), 18 deletions(-)

diff --git a/modules/uncertainty/R/compute_sobol_indices.R b/modules/uncertainty/R/compute_sobol_indices.R
index b6db1864712..154e9379c0f 100644
--- a/modules/uncertainty/R/compute_sobol_indices.R
+++ b/modules/uncertainty/R/compute_sobol_indices.R
@@ -7,28 +7,29 @@
 #' First-order indices quantify the share of output variance attributable to a
 #' parameter alone, while total-order indices summarize the full contribution of
 #' that parameter including its interactions with other parameters. PEcAn uses
-#' `sensobol` for these variance-based estimators; readers wanting estimator
-#' details or broader background should refer to Puy et al. (2022).
+#' `sensobol` for these variance-based estimators; see Saltelli et al. (2008)
+#' for methodological background and Puy et al. (2022) for package details.
 #'
 #' @param outdir PEcAn run output directory containing `ensemble.output.*.Rdata`
 #'   files.
 #' @param sobol_obj object produced by
 #'   `PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)`.
 #' @param var Variable name to summarize (default `"GPP"`).
-#' @param stat_fun Summary statistic applied to `var`. Retained for backwards
-#'   compatibility; standardized ensemble outputs are already scalar summaries.
 #'
 #' @return A tibble of Sobol first-order and total-order indices with attached
 #'   factor metadata.
-#' @references Puy, A., Lo Piano, S., Saltelli, A., and Levin, S. A. (2022).
+#' @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) {
+                                  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 ",
diff --git a/modules/uncertainty/man/compute_sobol_indices.Rd b/modules/uncertainty/man/compute_sobol_indices.Rd
index 4ebf55ece50..15fd8d350a8 100644
--- a/modules/uncertainty/man/compute_sobol_indices.Rd
+++ b/modules/uncertainty/man/compute_sobol_indices.Rd
@@ -4,23 +4,40 @@
 \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 `ensemble.output.*.Rdata`
+files.}
 
-\item{sobol_obj}{object produced by PEcAn.uncertainty::generate_joint_ensemble_design()}
+\item{sobol_obj}{object produced by
+`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 `"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 `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
+parameter alone, while total-order indices summarize the full contribution of
+that parameter including its interactions with other parameters. PEcAn uses
+`sensobol` for these variance-based estimators; see Saltelli et al. (2008)
+for methodological background and Puy et al. (2022) for package details.
+}
+\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}
 }

From 949edae2ed4ab2ab92c3e71eebbaf26583821cb4 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:43:23 -0700
Subject: [PATCH 11/25] Add error handling for input paths and validate sampled
 IDs in input.ens.gen function

---
 modules/uncertainty/R/ensemble.R | 27 ++++++++++++++++++++-------
 1 file changed, 20 insertions(+), 7 deletions(-)

diff --git a/modules/uncertainty/R/ensemble.R b/modules/uncertainty/R/ensemble.R
index 6863135ea32..841ce5d7a5a 100644
--- a/modules/uncertainty/R/ensemble.R
+++ b/modules/uncertainty/R/ensemble.R
@@ -562,14 +562,20 @@ input.ens.gen <- function(settings, ensemble_size, input, method = "sampling", p
 
   #-- assing the sample ids based on different scenarios
   input_path <- settings$run$inputs[[tolower(input)]]$path
+  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 - forexample, parent id may send id number 200 but we have only100 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))
+    # Sample replacement ids for inherited indices that exceed the child input bank.
+    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),
@@ -580,9 +586,16 @@ 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]
 
   return(samples)
 }
-

From 12e21a55a08d78ae2f11a0f12789e68b0a8c4ce0 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:44:28 -0700
Subject: [PATCH 12/25] ran scripts/generate_dependencies.R b/c change from
 sensitivity to sensobol package.

---
 docker/depends/pecan_package_dependencies.csv | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv
index d839a811865..3e00f6abb34 100644
--- a/docker/depends/pecan_package_dependencies.csv
+++ b/docker/depends/pecan_package_dependencies.csv
@@ -560,7 +560,7 @@
 "RPostgreSQL","*","base/db","Suggests",FALSE
 "RPostgreSQL","*","models/biocro","Suggests",FALSE
 "RSQLite","*","base/db","Suggests",FALSE
-"sensitivity","*","modules/uncertainty","Imports",FALSE
+"sensobol","*","modules/uncertainty","Imports",FALSE
 "sessioninfo","*","base/all","Suggests",FALSE
 "sf","*","modules/assim.sequential","Suggests",FALSE
 "sf","*","modules/data.atmosphere","Imports",FALSE

From 5f4ca3ed0797448525161a62e5a2805ed788ba43 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:47:46 -0700
Subject: [PATCH 13/25] Preserve parameter-parent input alignment in joint
 designs. Update documentation for generate_joint_ensemble_design: clarify
 parameters and return value

---
 .../R/generate_joint_ensemble_design.R        | 42 +++++++++++--------
 .../man/generate_joint_ensemble_design.Rd     | 18 ++++----
 2 files changed, 33 insertions(+), 27 deletions(-)

diff --git a/modules/uncertainty/R/generate_joint_ensemble_design.R b/modules/uncertainty/R/generate_joint_ensemble_design.R
index ed7da584791..82d97af6335 100644
--- a/modules/uncertainty/R/generate_joint_ensemble_design.R
+++ b/modules/uncertainty/R/generate_joint_ensemble_design.R
@@ -1,19 +1,3 @@
-#' 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.
-#'
-#' @param settings A PEcAn settings object containing ensemble configuration
-#' @param ensemble_size Integer specifying the number of ensemble members.
-#'   When `sobol = TRUE`, this is the Sobol base sample size `N`, not the
-#'   expanded number of model runs.
-#' @param sobol Logical, generate a variance-based Sobol design using
-#'   `sensobol`.
-#' @return A list with component `X`, a data frame design matrix describing
-#'   PEcAn parameter and sampled-input indices. If `sobol = TRUE`, the list
-#'   also includes the metadata needed by `compute_sobol_indices()`.
-#' @export
-
 .sobol_parameter_bank_size <- function(samples_file) {
   if (!file.exists(samples_file)) {
     return(0L)
@@ -44,6 +28,21 @@
   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.
+#'
+#' @param settings A PEcAn settings object containing ensemble configuration
+#' @param ensemble_size Integer specifying the number of ensemble members.
+#'   When `sobol = TRUE`, this is the Sobol base sample size `N`, not the
+#'   expanded number of model runs.
+#' @param sobol Logical, generate a variance-based Sobol design using
+#'   `sensobol`.
+#' @return A list with component `X`, a data frame design matrix describing
+#'   PEcAn parameter and sampled-input indices. If `sobol = TRUE`, the list
+#'   also includes the metadata needed by `compute_sobol_indices()`.
+#' @export
 generate_joint_ensemble_design <- function(settings,
                                            ensemble_size,
                                            sobol = FALSE) {
@@ -98,9 +97,13 @@ generate_joint_ensemble_design <- function(settings,
       sobol_design[["param"]],
       total_runs
     )
+    sampled_inputs[["parameters"]] <- list(ids = sobol_indices[["param"]])
 
     for (input_tag in setdiff(sobol_factors, "param")) {
       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)
@@ -154,8 +157,13 @@ generate_joint_ensemble_design <- function(settings,
   }
 
   # loop over inputs.
+  sampled_inputs[["parameters"]] <- list(ids = seq_len(ensemble_size))
   for (i in seq_along(samp.ordered)) {
     input_tag <- names(samp.ordered)[i]
+    if (identical(input_tag, "parameters")) {
+      next
+    }
+
     parent_name <- samp.ordered[[i]]$parent
 
     parent_ids <- if (!is.null(parent_name)) {
@@ -185,7 +193,7 @@ generate_joint_ensemble_design <- function(settings,
     )
   }
 
-  design_list[["param"]] <- seq_len(ensemble_size)
+  design_list[["param"]] <- sampled_inputs[["parameters"]]$ids
   design_matrix <- data.frame(design_list)
   return(list(X = design_matrix))
 }
diff --git a/modules/uncertainty/man/generate_joint_ensemble_design.Rd b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
index 717c96103fe..a0e6e446cee 100644
--- a/modules/uncertainty/man/generate_joint_ensemble_design.Rd
+++ b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
@@ -12,19 +12,17 @@ generate_joint_ensemble_design(settings, ensemble_size, sobol = FALSE)
 \arguments{
 \item{settings}{A PEcAn settings object containing ensemble configuration}
 
-\item{ensemble_size}{Integer specifying the number of ensemble members
-Since the `input_design` will only be generated once for the entire model run,
-the only situation, where we might want to recycle the existing `ensemble_samples`,
-is when we split and submit the larger SDA runs (e.g., 8,000 sites) into 
-smaller SDA experiments (e.g., 100 sites per job), where we want to keep using 
-the same parameters rather than creating new parameters for each job.}
+\item{ensemble_size}{Integer specifying the number of ensemble members.
+When `sobol = TRUE`, this is the Sobol base sample size `N`, not the
+expanded number of model runs.}
 
-\item{sobol}{for activating sobol}
+\item{sobol}{Logical, generate a variance-based Sobol design using
+`sensobol`.}
 }
 \value{
-A list containing ensemble samples and indices
-  If `sobol = TRUE`, the list will be a `sensitivity::soboljansen()` 
-  result and will contain the components documented therein.
+A list with component `X`, a data frame design matrix describing
+  PEcAn parameter and sampled-input indices. If `sobol = TRUE`, the list
+  also includes the metadata needed by `compute_sobol_indices()`.
 }
 \description{
 Generate joint ensemble design for parameter sampling

From 508b90a94be6f2e57a7898cff65b6e22048a4c93 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:52:02 -0700
Subject: [PATCH 14/25] Add Sobol tests for parameter-parent design behavior

---
 .../uncertainty/tests/testthat/test.sobol.R   | 64 +++++++++++++++++++
 1 file changed, 64 insertions(+)

diff --git a/modules/uncertainty/tests/testthat/test.sobol.R b/modules/uncertainty/tests/testthat/test.sobol.R
index e0e6af35ca2..6fcbddafd46 100644
--- a/modules/uncertainty/tests/testthat/test.sobol.R
+++ b/modules/uncertainty/tests/testthat/test.sobol.R
@@ -26,6 +26,30 @@ make_sobol_settings <- function(outdir) {
   )
 }
 
+make_parameter_parent_settings <- function(outdir, poolinit_paths) {
+  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(
+      poolinitcond = list(path = poolinit_paths)
+    )),
+    ensemble = list(
+      samplingspace = list(
+        parameters = list(method = "uniform"),
+        poolinitcond = list(method = "looping", parent = "parameters")
+      )
+    )
+  )
+}
+
 test_that("Sobol design expands to N * (k + 2) rows with metadata", {
   withr::with_tempdir({
     settings <- make_sobol_settings(getwd())
@@ -68,6 +92,46 @@ test_that("Non-Sobol design generation remains row-for-row", {
   })
 })
 
+test_that("Parameter-parented inputs inherit parameter ids in ordinary designs", {
+  withr::with_tempdir({
+    settings <- make_parameter_parent_settings(
+      getwd(),
+      paste0("ic", seq_len(8))
+    )
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 5,
+      sobol = FALSE
+    )
+
+    expect_equal(result$X$param, seq_len(5))
+    expect_equal(result$X$poolinitcond, result$X$param)
+  })
+})
+
+test_that("Sobol keeps parameter-parented child inputs in bounds", {
+  withr::with_tempdir({
+    set.seed(1)
+    settings <- make_parameter_parent_settings(
+      getwd(),
+      c("ic1", "ic2", "ic3")
+    )
+
+    result <- generate_joint_ensemble_design(
+      settings = settings,
+      ensemble_size = 4,
+      sobol = TRUE
+    )
+
+    expect_equal(nrow(result$X), 12)
+    expect_true(all(result$X$param >= 1))
+    expect_true(all(result$X$param <= 12))
+    expect_true(all(result$X$poolinitcond >= 1))
+    expect_true(all(result$X$poolinitcond <= 3))
+  })
+})
+
 test_that("compute_sobol_indices matches direct sensobol results", {
   withr::with_tempdir({
     sobol_obj <- list(

From e69e4f28c8d6b913374f0720fd2a99160fe2846e Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 17:55:04 -0700
Subject: [PATCH 15/25] Update uncertainty analysis tutorial to include global
 sensitivity analysis and clarify analysis types

---
 .../uncertainty.qmd                           | 113 ++++++++++++++----
 1 file changed, 92 insertions(+), 21 deletions(-)

diff --git a/documentation/tutorials/Demo_02_Uncertainty_Analysis/uncertainty.qmd b/documentation/tutorials/Demo_02_Uncertainty_Analysis/uncertainty.qmd
index caf15e3495e..e411109b6dc 100644
--- a/documentation/tutorials/Demo_02_Uncertainty_Analysis/uncertainty.qmd
+++ b/documentation/tutorials/Demo_02_Uncertainty_Analysis/uncertainty.qmd
@@ -14,13 +14,16 @@ format:
 
 # Introduction {#introduction}
 
-In Demo 2 we will be looking at how PEcAn can use information about parameter uncertainty to perform three automated analyses:
+In Demo 2 we will be looking at how PEcAn can use information about parameter uncertainty to perform four related analyses:
 
-- **Ensemble Analysis**: Repeat numerous model runs, each sampling from the parameter uncertainty, to generate a probability distribution of model projections. Allows us to put a confidence interval on the model.
-- **Sensitivity Analysis**: Repeats numerous model runs to assess how changes in model parameters will affect model outputs. Allows us to identify which parameters the model is most sensitive to.
-- **Uncertainty Analysis**: Combines information about model sensitivity with information about parameter uncertainty to determine the contribution of each model parameter to the uncertainty in model outputs. Allow us to identify which parameters are driving model uncertainty.
+- **Ensemble Analysis**: Repeats model runs while sampling from parameter uncertainty to propagate parameter uncertainty through to model output and estimate uncertainty in model projections (LeBauer et al. 2013).
+- **Sensitivity Analysis**: Repeats model runs to assess how changes in model parameters affect model outputs. Allows us to identify which parameters the model is most sensitive to.
+- **Variance Decomposition**: Combines the one-at-a-time sensitivity results with parameter uncertainty to estimate which parameters contribute most to model uncertainty and help prioritize subsequent data collection (LeBauer et al. 2013).
+- **Global Sensitivity Analysis**: Uses a Sobol design to partition output variance across major sources of uncertainty such as parameters, meteorology, initial conditions, and events. In contrast to one-at-a-time analysis, this global design can account for interactions among factors (Saltelli et al. 2008).
 
-This demo shows how to run an uncertainty analysis workflow in PEcAn using an R-based Quarto notebook. It covers loading settings, configuring models, running simulations, and performing ensemble and sensitivity analyses to assess uncertainty and parameter importance. This programmatic approach complements the web-based PEcAn interface.
+This demo shows how to run an uncertainty analysis workflow in PEcAn using an R-based Quarto notebook. It covers loading settings, configuring models, running simulations, and performing ensemble analysis, one-at-a-time sensitivity analysis, variance decomposition, and Sobol global sensitivity analysis. This programmatic approach complements the web-based PEcAn interface.
+
+Dietze (2017) provides a comprehensive review of sensitivity and uncertainty analysis in ecosystem modeling.
 
 **Context & modeling scenario:**
 
@@ -30,7 +33,7 @@ We simulate plant and ecosystem carbon balance (Net Primary Productivity and Net
 
 1. Configure a PEcAn workflow by loading and validating a `pecan.xml` settings file.
 2. Run a set of ecosystem model simulations by writing model configuration files and then running SIPNET.
-3. Quantify uncertainty using ensemble analysis, sensitivity analyses, and variance decomposition.
+3. Quantify uncertainty using ensemble analysis, one-at-a-time sensitivity analysis, variance decomposition, and Sobol global sensitivity analysis.
 4. Visualize results to identify important parameters and how they influence model variance.
 5. Change configuration settings and re-run the workflow.
 
@@ -102,7 +105,7 @@ The sampling method for generating parameter sets is set to **halton** in this d
 
 PEcAn's sensitivity analysis includes a handy shortcut that converts a specified standard deviation into its normal quantile equivalent. In the example pecan.xml, these are set to **-1, 1** (the median value, 0, occurs by default) which are converted internally to the 15.9th and 84.1th quantiles of the parameter distribution. You can add more quantiles to explore a wider range of parameter values: {-3, -2, -1, 1, 2, 3} is often used in practice. 
 
-By working in quantiles relative to each parameter’s distribution, the sensitivity and variance decomposition reflect sensitivity across the range of parameter values. Many sensitivity analyses tools use a fixed perturbation size such as the mean +/- 10%. PEcAn's SA does not take this approach because it does not capture the uncertainty across the parameter distribution and can not be used for variance decomposition.
+By working in quantiles relative to each parameter's distribution, the sensitivity and variance decomposition reflect sensitivity across the range of parameter values. This matches the PEcAn goal of propagating parameter uncertainty through to model output and then evaluating how much each parameter contributes to model uncertainty (LeBauer et al. 2013). Many one-at-a-time analyses use a fixed perturbation such as mean +/- 10%, but that does not capture uncertainty across the parameter distribution and therefore does not support the same variance decomposition workflow.
 
 ## Load the settings file
 
@@ -149,9 +152,9 @@ Next we convert all of the model output from the previous run to a standard form
 runModule.get.results(settings)
 ```
 
-# Ensemble and Sensitivity Analysis {#sec-run-uncertainty}
+# Ensemble, Sensitivity, Variance Decomposition, and Global Sensitivity Analysis {#sec-run-uncertainty}
 
-Next, we use the outputs from the previous step to perform ensemble and sensitivity analyses.
+Next, we use the outputs from the previous step for two related workflows: an ensemble plus one-at-a-time sensitivity workflow that feeds the variance decomposition analysis, and a separate Sobol global sensitivity workflow.
 
 **Ensemble Analysis**: Quantifies uncertainty in model predictions by running multiple simulations with parameters sampled from their uncertainty distributions. This generates probability distributions of model outputs and confidence intervals.
 
@@ -159,19 +162,54 @@ Next, we use the outputs from the previous step to perform ensemble and sensitiv
 runModule.run.ensemble.analysis(settings)
 ```
 
-**Sensitivity Analysis**: Systematically varies individual parameters to assess their influence on model outputs. This identifies which parameters most strongly affect model predictions and helps prioritize parameter refinement efforts.
+**Sensitivity Analysis**: Systematically varies individual parameters to assess their influence on model outputs. In this demo, these one-at-a-time runs are also the basis for the variance decomposition analysis that follows. This is useful for prioritizing parameters, but it does not estimate interactions among uncertain factors.
 
 ```{r run-sensitivity-analysis}
 runModule.run.sensitivity.analysis(settings)
 ```
 
+## Sobol Global Sensitivity Analysis
+
+Sobol global sensitivity analysis is a separate variance-based workflow. Rather than using the one-at-a-time sensitivity results, it starts from a Sobol design and estimates first-order and total-order contributions directly from a global run ensemble. This is a variance-based global sensitivity analysis: it diagnoses how uncertainty in the inputs contributes to uncertainty in the model output, including interactions among factors (Dietze 2017; Saltelli et al. 2008). PEcAn computes these indices using `sensobol` (Puy et al. 2022). It is often helpful to run this analysis in a new output directory to clearly separate runs generated by the Sobol design.
+
+> Note: In Sobol mode, `ensemble_size` means the base sample size `N`, not the expanded number of model runs.
+
+```{r run-sobol-gsa, eval=FALSE}
+gsa_settings <- settings
+# Use a fresh output directory before running this block.
+
+sobol_design <- PEcAn.uncertainty::generate_joint_ensemble_design(
+  settings = gsa_settings,
+  ensemble_size = 16,
+  sobol = TRUE
+)
+
+gsa_settings <- PEcAn.workflow::runModule.run.write.configs(
+  settings = gsa_settings,
+  input_design = sobol_design$X
+)
+PEcAn.workflow::runModule_start_model_runs(gsa_settings)
+runModule.get.results(gsa_settings)
+
+sobol_indices <- PEcAn.uncertainty::compute_sobol_indices(
+  outdir = gsa_settings$outdir,
+  sobol_obj = sobol_design,
+  var = "NPP"
+)
+
+sobol_indices
+```
+
 # PEcAn Outputs {#sec-outputs}
 
 ## Output Directory Structure
 
 These are the key folders and files that will be created under the directory defined by `settings$outdir` (e.g., `demo_outdir` in the example). The file contents are described in the next section.
 
-We discussed the output directory in Demo 1 (Basic Run), but now we have three new folders that contain outputs from the sensitivity, ensemble, and variance decomposition analyses. Here we focus on the additional outputs generated by the ensemble and sensitivity analyses.
+We discussed the output directory in Demo 1 (Basic Run), but now we have additional outputs from the ensemble, one-at-a-time sensitivity, variance decomposition, and Sobol global sensitivity analyses. Here we focus on those additional analysis products.
+
+Note: The layout below is representative - the actual files in the out directory will depend
+on workflow settings.
 
 ```
 demo_outdir/
@@ -186,6 +224,8 @@ demo_outdir/
 ├── ensemble.samples.*.Rdata              # Parameter samples used for ensemble
 ├── ensemble.ts.*.Rdata                   # Time series data from ensemble
 ├── samples.Rdata                         # Parameter samples for both SA and ensemble
+├── sobol.design.*.Rdata                  # Saved Sobol design metadata for post-processing
+├── sobol.indices.*.Rdata                 # Saved Sobol first-order and total-order indices
 ├── sensitivity.output.*.Rdata            # SA model outputs
 ├── sensitivity.results.*.Rdata           # Processed SA results
 ├── sensitivity.samples.*.Rdata           # SA parameter samples
@@ -204,7 +244,7 @@ demo_outdir/
 
 # Understanding PEcAn Uncertainty Analysis Outputs
 
-After running ensemble and sensitivity analyses, PEcAn generates several important outputs that help you understand model uncertainty and parameter sensitivity.
+After running ensemble analysis, one-at-a-time sensitivity analysis, variance decomposition, and global sensitivity analysis, PEcAn generates several important outputs that help you understand model uncertainty and parameter sensitivity.
 
 The `samples.Rdata` file contains the parameter values used in the sensitivity and ensemble runs. It stores two objects, `sa.samples` and `ensemble.samples`, which are the parameter values for the sensitivity analysis and ensemble runs, respectively.
 
@@ -212,7 +252,7 @@ The `samples.Rdata` file contains the parameter values used in the sensitivity a
 
 The ensemble analysis produces:
 
-- **`ensemble.Rdata`**: Contains `ensemble.output` object with model predictions for all ensemble members
+- **`ensemble.output.[RunID].[Variable].[StartYear].[EndYear].Rdata`**: Contains the `ensemble.output` object with model predictions for all ensemble members
 - **`ensemble.analysis.[RunID].[Variable].[StartYear].[EndYear].pdf`**: Histogram and boxplot of ensemble predictions
 - **`ensemble.ts.[RunID].[Variable].[StartYear].[EndYear].pdf`**: Time-series plot showing ensemble mean, median, and 95% confidence intervals
 
@@ -233,6 +273,16 @@ The variance decomposition produces:
   - Elasticity (normalized sensitivity) 
   - Partial standard deviation of each parameter
 
+
+## Sobol Global Sensitivity Outputs
+
+The Sobol global sensitivity analysis generates:
+
+- **`sobol.design.[RunID].[Variable].[StartYear].[EndYear].Rdata`**: The saved Sobol design metadata used to map PEcAn runs back to the Sobol estimator.
+- **`sobol.indices.[RunID].[Variable].[StartYear].[EndYear].Rdata`**: The computed first-order and total-order Sobol indices for each uncertainty class in the design. First-order indices summarize the direct contribution of each factor, whereas total-order indices summarize the factor's overall contribution including interactions (Saltelli et al. 2008). PEcAn computes these indices using `sensobol` (Puy et al. 2022).
+
+After running the Sobol workflow, you can inspect `sobol_indices` directly or load the saved `sobol.indices.*.Rdata` file to compare first-order and total-order contributions across uncertainty classes.
+
 ## Interpreting the Results
 
 **Variance Decomposition Analysis:**
@@ -268,7 +318,7 @@ variable <- settings$ensemble$variable
 pft <- settings$pfts[[1]]
 start.year <- lubridate::year(settings$run$start.date)
 end.year <- lubridate::year(settings$run$end.date)
-ensemble.id <- if (!is.null(settings$ensemble$id)) settings$ensemble$id else "NOENSEMBLEID"
+ensemble.id <- if (!is.null(settings$ensemble$ensemble.id)) settings$ensemble$ensemble.id else "NOENSEMBLEID"
 
 # --- 2. Load and Plot Ensemble Output Distribution ---
 # This section visualizes the distribution of the ensemble model runs.
@@ -325,9 +375,14 @@ This block visualizes the results of the sensitivity analysis. The plots show ho
 # outputs needed to generate the sensitivity and variance decomposition plots.
 
 # Construct the path to the sensitivity analysis results file
+sensitivity.id <- if (!is.null(settings$sensitivity.analysis$ensemble.id)) {
+  settings$sensitivity.analysis$ensemble.id
+} else {
+  "NOENSEMBLEID"
+}
 sens_file <- file.path(
   settings$outdir,
-  paste0("sensitivity.results", ".", ensemble.id, ".", variable, ".", start.year, ".", end.year, ".Rdata")
+  paste0("sensitivity.results", ".", sensitivity.id, ".", variable, ".", start.year, ".", end.year, ".Rdata")
 )
 
 # Check if the file exists, then load the data and generate plots
@@ -453,7 +508,7 @@ plot(model_output$posix, model_output$GPP,
   col = "green",
   xlab = "Date",
   ylab = "Carbon Flux (kg C m-2 s-1)",
-  main = paste("Carbon Fluxes Over Time — PEcAn", runid)
+  main = paste("Carbon Fluxes Over Time - PEcAn", runid)
 )
 lines(model_output$posix, model_output$NPP, col = "blue")
 legend("topright", legend = c("GPP", "NPP"), col = c("green", "blue"), lty = 1)
@@ -468,7 +523,7 @@ plot(model_output$posix, model_output$TotLivBiom,
   col = "darkgreen",
   xlab = "Date",
   ylab = "Carbon Pool (kg C m-2)",
-  main = paste("Carbon Pools Over Time — PEcAn", runid)
+  main = paste("Carbon Pools Over Time - PEcAn", runid)
 )
 lines(model_output$posix, model_output$TotSoilCarb, col = "brown")
 legend("topright", legend = c("Total Live Biomass", "Total Soil Carbon"), col = c("darkgreen", "brown"), lty = 1)
@@ -483,7 +538,7 @@ plot(model_output$posix, model_output$SoilMoist,
   col = "blue",
   xlab = "Date",
   ylab = "Soil Moisture (kg m-2)",
-  main = paste("Soil Moisture Over Time — PEcAn", runid)
+  main = paste("Soil Moisture Over Time - PEcAn", runid)
 )
 lines(model_output$posix, model_output$SWE, col = "lightblue")
 legend("topright", legend = c("Soil Moisture", "Snow Water Equivalent"), col = c("blue", "lightblue"), lty = 1)
@@ -498,18 +553,24 @@ plot(model_output$posix, model_output$LAI,
   col = "darkgreen",
   xlab = "Date",
   ylab = "LAI (m2 m-2)",
-  main = paste("Leaf Area Index Over Time — PEcAn", runid)
+  main = paste("Leaf Area Index Over Time - PEcAn", runid)
 )
 lines(model_output$posix, model_output$AbvGrndWood, col = "brown")
 legend("topright", legend = c("LAI", "Above Ground Wood"), col = c("darkgreen", "brown"), lty = 1)
 ```
 
-# Conclusion
+# Your turn
 
-This notebook demonstrated how to set up, run, and analyze a PEcAn ecosystem model workflow programmatically. You can now modify parameters, try different models, or extend the analysis as needed.
+Try visualizing other variables in the model output! Use the `available_vars` object to see what variables are available and create your own time series plots. You can also customize the plots with different colors, labels, and titles.
 
 Try editing the `pecan.xml` file. Give it a new name and update the **settings_path** variable at the beginning of this Demo to point to the new file. See how the changes affect the model output!
 
+Inspect `sobol_indices` and identify which uncertainty class has the largest first-order contribution to variance, and which has the largest total-order contribution. Differences between first-order and total-order indices indicate the importance of interactions. <!-- TODO: add Sobol GSA visualization and interpretation example -->
+
+# Conclusion
+
+This notebook demonstrated how to set up, run, and analyze a PEcAn ecosystem model workflow programmatically. You can now modify parameters, try different models, or extend the analysis as needed.
+
 # Further Exploration
 
 The [next set of tutorials](#demo-table) will focus on the process of data assimilation and parameter estimation. The next two steps are in “.Rmd” files which can be viewed online. 
@@ -553,3 +614,13 @@ PEcAn.all::pecan_version()
 ```{r session-info}
 sessionInfo()
 ```
+
+# References
+
+Dietze, M., 2017. Ecological Forecasting. Princeton University Press.
+
+LeBauer, D. S., Wang, D., Richter, K. T., Davidson, C. C., and Dietze, M. C. 2013. Facilitating feedbacks between field measurements and ecosystem models. *Ecological Monographs* 83(2):133-154.
+
+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. https://doi.org/10.18637/jss.v102.i05
+
+Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D. et al. 2008. *Global Sensitivity Analysis: The Primer*. New York: John Wiley & Sons.

From 60cf71e41977c325c3ed4399c382445c87e318e6 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 19:14:53 -0700
Subject: [PATCH 16/25] Add sensobol changelog notes

---
 CHANGELOG.md                | 4 ++++
 modules/uncertainty/NEWS.md | 2 ++
 2 files changed, 6 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9b14f4ebe4e..2f770868d06 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).
@@ -18,6 +19,9 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Added datasets to `PEcAn.data.land` 
   *  `landiq_crop_mapping_codes` dataset mapping LandIQ crop classification codes to human-readable crop names.
   *  `bism_kc_by_crop` dataset containing BISm crop coefficient schedules and stage timing references for use in ET estimation, including columns that map to LandIQ class and subclass.
+
+### Fixed
+- Fixed joint ensemble design and workflow configuration handling to preserve parent-child input alignment and validate `input_design$param` indices before writing configs.
 - PEcAn.SIPNET gains support for SIPNET v2, whose features includes management events, nitrogen cycle tracking, explicit N2O and methane fluxes, runtime setting of feature flags, and changes to the parameter set (now 73 parameters). SIPNET v1 is still fully supported, but workarounds for bugs in the legacy `sipnet.unk` version have been removed.
 
 ### Fixed
diff --git a/modules/uncertainty/NEWS.md b/modules/uncertainty/NEWS.md
index 5a446773800..92cab17f25c 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.
+* Fixed parent-linked input design handling and validation so parameter-selected sample indices stay aligned and in bounds during workflow configuration.
 * 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

From 8dbf9beffd8dee07fbd330aa04ccf908f3534971 Mon Sep 17 00:00:00 2001
From: David LeBauer <dlebauer@gmail.com>
Date: Wed, 18 Mar 2026 19:34:53 -0700
Subject: [PATCH 17/25] Add trait sample bank size validation and related tests

---
 base/workflow/R/run.write.configs.R           | 44 ++++++++++++-
 .../testthat/test.run.write.configs.sobol.R   | 49 +++++++++++++++
 models/sipnet/DESCRIPTION                     |  1 -
 modules/uncertainty/DESCRIPTION               |  3 +-
 .../R/generate_joint_ensemble_design.R        | 44 ++++++++++---
 .../uncertainty/tests/testthat/test.sobol.R   | 61 +++++++++++++++++++
 6 files changed, 187 insertions(+), 15 deletions(-)

diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R
index 0ef132321de..144120f427c 100644
--- a/base/workflow/R/run.write.configs.R
+++ b/base/workflow/R/run.write.configs.R
@@ -27,6 +27,40 @@
 #'
 #' @author David LeBauer, Shawn Serbin, Ryan Kelly, Mike Dietze, Akash B V
 
+.trait_sample_bank_size <- function(trait.samples) {
+  if (is.null(trait.samples) || length(trait.samples) == 0) {
+    return(0L)
+  }
+
+  bank_sizes <- unlist(
+    lapply(trait.samples, function(pft_traits) {
+      if (is.null(pft_traits) || length(pft_traits) == 0) {
+        return(integer(0))
+      }
+
+      vapply(
+        pft_traits,
+        function(trait_values) {
+          if (is.null(trait_values) || length(trait_values) == 0) {
+            return(NA_integer_)
+          }
+
+          as.integer(length(trait_values))
+        },
+        integer(1)
+      )
+    }),
+    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))
+}
+
 run.write.configs <- function(settings, ensemble.size, input_design, write = TRUE,
                               posterior.files = rep(NA, length(settings$pfts)),
                               overwrite = TRUE) {
@@ -142,13 +176,17 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
           "samples.Rdata does not contain trait.samples required for input_design$param"
         )
       }
-      first_pft <- trait.samples[[1]]
-      first_trait <- first_pft[[1]]
       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")
       }
-      if (is.null(first_trait) || any(trait_sample_indices > length(first_trait))) {
+      parameter_bank_size <- .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"
         )
diff --git a/base/workflow/tests/testthat/test.run.write.configs.sobol.R b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
index 3320b29a0a9..c66b2173a25 100644
--- a/base/workflow/tests/testthat/test.run.write.configs.sobol.R
+++ b/base/workflow/tests/testthat/test.run.write.configs.sobol.R
@@ -32,3 +32,52 @@ test_that("runModule.run.write.configs uses input_design row count", {
   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/models/sipnet/DESCRIPTION b/models/sipnet/DESCRIPTION
index 2cc2896bc25..6936e013d59 100644
--- a/models/sipnet/DESCRIPTION
+++ b/models/sipnet/DESCRIPTION
@@ -14,7 +14,6 @@ URL: https://pecanproject.github.io
 BugReports: https://github.com/PecanProject/pecan/issues
 Depends: R (>= 4.1.0)
 Imports:
-    data.table,
     dplyr,
     jsonlite,
     lubridate (>= 1.6.0),
diff --git a/modules/uncertainty/DESCRIPTION b/modules/uncertainty/DESCRIPTION
index a7a5d34d216..d7973645550 100644
--- a/modules/uncertainty/DESCRIPTION
+++ b/modules/uncertainty/DESCRIPTION
@@ -43,7 +43,8 @@ Imports:
     purrr,
     randtoolbox,
     rlang,
-    sensobol
+    sensobol,
+    tibble
 Suggests:
     mockery,
     testthat (>= 1.0.2),
diff --git a/modules/uncertainty/R/generate_joint_ensemble_design.R b/modules/uncertainty/R/generate_joint_ensemble_design.R
index 6a46738b3bd..c412e131bf2 100644
--- a/modules/uncertainty/R/generate_joint_ensemble_design.R
+++ b/modules/uncertainty/R/generate_joint_ensemble_design.R
@@ -1,26 +1,50 @@
-.sobol_parameter_bank_size <- function(samples_file) {
-  if (!file.exists(samples_file)) {
+.trait_sample_bank_size <- function(trait.samples) {
+  if (is.null(trait.samples) || length(trait.samples) == 0) {
     return(0L)
   }
 
-  samples <- new.env(parent = emptyenv())
-  load(samples_file, envir = samples)
+  bank_sizes <- unlist(
+    lapply(trait.samples, function(pft_traits) {
+      if (is.null(pft_traits) || length(pft_traits) == 0) {
+        return(integer(0))
+      }
 
-  if (is.null(samples$trait.samples) || length(samples$trait.samples) == 0) {
+      vapply(
+        pft_traits,
+        function(trait_values) {
+          if (is.null(trait_values) || length(trait_values) == 0) {
+            return(NA_integer_)
+          }
+
+          as.integer(length(trait_values))
+        },
+        integer(1)
+      )
+    }),
+    use.names = FALSE
+  )
+
+  bank_sizes <- bank_sizes[!is.na(bank_sizes) & bank_sizes > 0L]
+  if (length(bank_sizes) == 0) {
     return(0L)
   }
 
-  first_pft <- samples$trait.samples[[1]]
-  if (is.null(first_pft) || length(first_pft) == 0) {
+  as.integer(min(bank_sizes))
+}
+
+.sobol_parameter_bank_size <- function(samples_file) {
+  if (!file.exists(samples_file)) {
     return(0L)
   }
 
-  first_trait <- first_pft[[1]]
-  if (is.null(first_trait)) {
+  samples <- new.env(parent = emptyenv())
+  load(samples_file, envir = samples)
+
+  if (is.null(samples$trait.samples) || length(samples$trait.samples) == 0) {
     return(0L)
   }
 
-  as.integer(length(first_trait))
+  .trait_sample_bank_size(samples$trait.samples)
 }
 
 .map_sobol_to_indices <- function(x, size) {
diff --git a/modules/uncertainty/tests/testthat/test.sobol.R b/modules/uncertainty/tests/testthat/test.sobol.R
index 6fcbddafd46..1f50a80e31b 100644
--- a/modules/uncertainty/tests/testthat/test.sobol.R
+++ b/modules/uncertainty/tests/testthat/test.sobol.R
@@ -50,6 +50,33 @@ make_parameter_parent_settings <- function(outdir, poolinit_paths) {
   )
 }
 
+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 expands to N * (k + 2) rows with metadata", {
   withr::with_tempdir({
     settings <- make_sobol_settings(getwd())
@@ -132,6 +159,40 @@ test_that("Sobol keeps parameter-parented child inputs in bounds", {
   })
 })
 
+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
+    )
+
+    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(

From 437a34f6827a6c490f88bdff545ee03d2b484d95 Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:15:27 -0400
Subject: [PATCH 18/25] make events an independent Sobol factor, remove
 parent-child filter

---
 modules/uncertainty/R/compute_sobol_indices.R |  30 +++-
 .../R/generate_joint_ensemble_design.R        | 135 +++++++-----------
 .../uncertainty/tests/testthat/test.sobol.R   | 111 +++++---------
 3 files changed, 114 insertions(+), 162 deletions(-)

diff --git a/modules/uncertainty/R/compute_sobol_indices.R b/modules/uncertainty/R/compute_sobol_indices.R
index 154e9379c0f..e30e62459c0 100644
--- a/modules/uncertainty/R/compute_sobol_indices.R
+++ b/modules/uncertainty/R/compute_sobol_indices.R
@@ -1,20 +1,36 @@
 #' Compute Sobol indices from a finished PEcAn run
 #'
 #' Loads standardized ensemble output for a Sobol run, computes first-order and
-#' total-order Sobol indices with `sensobol`, and saves both the Sobol design
+#' total-order Sobol indices with \code{sensobol}, and saves both the Sobol design
 #' metadata and the computed indices using PEcAn-style filenames.
 #'
 #' First-order indices quantify the share of output variance attributable to a
-#' parameter alone, while total-order indices summarize the full contribution of
-#' that parameter including its interactions with other parameters. PEcAn uses
-#' `sensobol` for these variance-based estimators; see Saltelli et al. (2008)
+#' 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.
 #'
-#' @param outdir PEcAn run output directory containing `ensemble.output.*.Rdata`
+#' 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
-#'   `PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)`.
-#' @param var Variable name to summarize (default `"GPP"`).
+#'   \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.
diff --git a/modules/uncertainty/R/generate_joint_ensemble_design.R b/modules/uncertainty/R/generate_joint_ensemble_design.R
index c412e131bf2..5967ee1ce8a 100644
--- a/modules/uncertainty/R/generate_joint_ensemble_design.R
+++ b/modules/uncertainty/R/generate_joint_ensemble_design.R
@@ -1,25 +1,30 @@
-.trait_sample_bank_size <- function(trait.samples) {
+#' 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(
-    lapply(trait.samples, function(pft_traits) {
+    purrr::map(trait.samples, function(pft_traits) {
       if (is.null(pft_traits) || length(pft_traits) == 0) {
         return(integer(0))
       }
-
-      vapply(
-        pft_traits,
-        function(trait_values) {
-          if (is.null(trait_values) || length(trait_values) == 0) {
-            return(NA_integer_)
-          }
-
-          as.integer(length(trait_values))
-        },
-        integer(1)
-      )
+      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
   )
@@ -44,7 +49,7 @@
     return(0L)
   }
 
-  .trait_sample_bank_size(samples$trait.samples)
+  trait_sample_bank_size(samples$trait.samples)
 }
 
 .map_sobol_to_indices <- function(x, size) {
@@ -53,21 +58,31 @@
 }
 
 #' 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.
 #'
+#' 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.
+#'
 #' @param settings A PEcAn settings object containing ensemble configuration.
 #' @param ensemble_size Integer specifying the number of ensemble members.
 #'   When \code{sobol = TRUE}, this is the Sobol base sample size \code{N}, not
-#'   the expanded number of model runs.
+#'   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 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{\link{compute_sobol_indices}}: \code{N}, \code{params},
+#'   \code{backend}, \code{matrices}, \code{first}, \code{total}, and
+#'   \code{factor_metadata}.
 #'
 #' @export
 
@@ -77,28 +92,14 @@ generate_joint_ensemble_design <- function(settings,
   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)])]
 
   if (sobol) {
-    sobol_factors <- c(
-      "param",
-      names(samp.ordered)[
-        names(samp.ordered) != "parameters" &
-          vapply(
-            samp.ordered,
-            function(x) is.null(x$parent),
-            logical(1)
-          )
-      ]
-    )
+    # 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")
@@ -120,6 +121,7 @@ generate_joint_ensemble_design <- function(settings,
     )
     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"]],
@@ -127,49 +129,31 @@ generate_joint_ensemble_design <- function(settings,
     )
     sampled_inputs[["parameters"]] <- list(ids = sobol_indices[["param"]])
 
-    for (input_tag in setdiff(sobol_factors, "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")
+        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)
       )
-    }
-
-    for (i in seq_along(samp.ordered)) {
-      input_tag <- names(samp.ordered)[i]
-
-      if (identical(input_tag, "parameters")) {
-        next
-      }
-
-      parent_name <- samp.ordered[[i]]$parent
-      if (!is.null(parent_name)) {
-        input_result <- PEcAn.uncertainty::input.ens.gen(
-          settings = settings,
-          ensemble_size = total_runs,
-          input = input_tag,
-          method = samp.ordered[[i]]$method,
-          parent_ids = sampled_inputs[[parent_name]]
-        )
-        sampled_inputs[[input_tag]] <- input_result
-        design_list[[input_tag]] <- input_result$ids
-      } else if (input_tag %in% names(sobol_indices)) {
-        sampled_inputs[[input_tag]] <- list(ids = sobol_indices[[input_tag]])
-        design_list[[input_tag]] <- sobol_indices[[input_tag]]
-      }
+      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 <- data.frame(design_list)
+    design_matrix <- tibble::as_tibble(design_list)
 
-    factor_metadata <- data.frame(
+    factor_metadata <- tibble::tibble(
       factor = sobol_factors,
       source_type = sobol_factors,
-      source_tag = ifelse(sobol_factors == "param", NA_character_, sobol_factors),
-      stringsAsFactors = FALSE
+      source_tag = ifelse(
+        sobol_factors == "param", NA_character_, sobol_factors
+      )
     )
 
     return(list(
@@ -184,28 +168,15 @@ generate_joint_ensemble_design <- function(settings,
     ))
   }
 
+  # non-Sobol path: simple sequential or sampled design
   sampled_inputs[["parameters"]] <- list(ids = seq_len(ensemble_size))
-  for (i in seq_along(samp.ordered)) {
-    input_tag <- names(samp.ordered)[i]
-    if (identical(input_tag, "parameters")) {
-      next
-    }
-
-    parent_name <- samp.ordered[[i]]$parent
-    parent_ids <- if (!is.null(parent_name)) {
-      sampled_inputs[[parent_name]]
-    } else {
-      NULL
-    }
-
+  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
     design_list[[input_tag]] <- input_result$ids
   }
@@ -220,6 +191,6 @@ generate_joint_ensemble_design <- function(settings,
   }
 
   design_list[["param"]] <- sampled_inputs[["parameters"]]$ids
-  design_matrix <- data.frame(design_list)
+  design_matrix <- tibble::as_tibble(design_list)
   return(list(X = design_matrix))
 }
diff --git a/modules/uncertainty/tests/testthat/test.sobol.R b/modules/uncertainty/tests/testthat/test.sobol.R
index 1f50a80e31b..4ffcea7aa42 100644
--- a/modules/uncertainty/tests/testthat/test.sobol.R
+++ b/modules/uncertainty/tests/testthat/test.sobol.R
@@ -20,31 +20,7 @@ make_sobol_settings <- function(outdir) {
         parameters = list(method = "uniform"),
         met = list(method = "sampling"),
         poolinitcond = list(method = "looping"),
-        events = list(method = "sampling", parent = "met")
-      )
-    )
-  )
-}
-
-make_parameter_parent_settings <- function(outdir, poolinit_paths) {
-  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(
-      poolinitcond = list(path = poolinit_paths)
-    )),
-    ensemble = list(
-      samplingspace = list(
-        parameters = list(method = "uniform"),
-        poolinitcond = list(method = "looping", parent = "parameters")
+        events = list(method = "sampling")
       )
     )
   )
@@ -77,7 +53,7 @@ make_mixed_bank_sobol_settings <- function(outdir) {
   )
 }
 
-test_that("Sobol design expands to N * (k + 2) rows with metadata", {
+test_that("Sobol design treats all inputs as independent factors", {
   withr::with_tempdir({
     settings <- make_sobol_settings(getwd())
 
@@ -87,19 +63,46 @@ test_that("Sobol design expands to N * (k + 2) rows with metadata", {
       sobol = TRUE
     )
 
-    expect_equal(nrow(result$X), 20)
+    # 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"))
+    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")
-    expect_true(all(c("param", "met", "poolinitcond", "events") %in% names(result$X)))
-    expect_equal(result$X$events, result$X$met)
+
+    # 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 <= 20))
-    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"))
+    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")
+    )
   })
 })
 
@@ -119,46 +122,6 @@ test_that("Non-Sobol design generation remains row-for-row", {
   })
 })
 
-test_that("Parameter-parented inputs inherit parameter ids in ordinary designs", {
-  withr::with_tempdir({
-    settings <- make_parameter_parent_settings(
-      getwd(),
-      paste0("ic", seq_len(8))
-    )
-
-    result <- generate_joint_ensemble_design(
-      settings = settings,
-      ensemble_size = 5,
-      sobol = FALSE
-    )
-
-    expect_equal(result$X$param, seq_len(5))
-    expect_equal(result$X$poolinitcond, result$X$param)
-  })
-})
-
-test_that("Sobol keeps parameter-parented child inputs in bounds", {
-  withr::with_tempdir({
-    set.seed(1)
-    settings <- make_parameter_parent_settings(
-      getwd(),
-      c("ic1", "ic2", "ic3")
-    )
-
-    result <- generate_joint_ensemble_design(
-      settings = settings,
-      ensemble_size = 4,
-      sobol = TRUE
-    )
-
-    expect_equal(nrow(result$X), 12)
-    expect_true(all(result$X$param >= 1))
-    expect_true(all(result$X$param <= 12))
-    expect_true(all(result$X$poolinitcond >= 1))
-    expect_true(all(result$X$poolinitcond <= 3))
-  })
-})
-
 test_that("Sobol regenerates parameter bank when any PFT bank is too short", {
   withr::with_tempdir({
     settings <- make_mixed_bank_sobol_settings(getwd())
@@ -181,6 +144,8 @@ test_that("Sobol regenerates parameter bank when any PFT bank is too short", {
       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,

From 7469668761bc154320e6ad60036d10038ef73b75 Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:16:18 -0400
Subject: [PATCH 19/25] remove redundent funciton signature and made it export

---
 base/workflow/R/run.write.configs.R | 36 +----------------------------
 1 file changed, 1 insertion(+), 35 deletions(-)

diff --git a/base/workflow/R/run.write.configs.R b/base/workflow/R/run.write.configs.R
index 144120f427c..4d39fd3c37e 100644
--- a/base/workflow/R/run.write.configs.R
+++ b/base/workflow/R/run.write.configs.R
@@ -27,40 +27,6 @@
 #'
 #' @author David LeBauer, Shawn Serbin, Ryan Kelly, Mike Dietze, Akash B V
 
-.trait_sample_bank_size <- function(trait.samples) {
-  if (is.null(trait.samples) || length(trait.samples) == 0) {
-    return(0L)
-  }
-
-  bank_sizes <- unlist(
-    lapply(trait.samples, function(pft_traits) {
-      if (is.null(pft_traits) || length(pft_traits) == 0) {
-        return(integer(0))
-      }
-
-      vapply(
-        pft_traits,
-        function(trait_values) {
-          if (is.null(trait_values) || length(trait_values) == 0) {
-            return(NA_integer_)
-          }
-
-          as.integer(length(trait_values))
-        },
-        integer(1)
-      )
-    }),
-    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))
-}
-
 run.write.configs <- function(settings, ensemble.size, input_design, write = TRUE,
                               posterior.files = rep(NA, length(settings$pfts)),
                               overwrite = TRUE) {
@@ -180,7 +146,7 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
       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 <- .trait_sample_bank_size(trait.samples)
+      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"

From 1ba2eefc99ff8c47b4d33b2d9c4159511224c4ef Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:17:16 -0400
Subject: [PATCH 20/25] update NAMESPACE

---
 modules/uncertainty/NAMESPACE | 1 +
 1 file changed, 1 insertion(+)

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,"%>%")

From 7afbb7ef157b558b5d41fce58ba78734ddcee696 Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:17:47 -0400
Subject: [PATCH 21/25] update compute_sobol_indices.Rd

---
 .../uncertainty/man/compute_sobol_indices.Rd  | 31 ++++++++++++++-----
 1 file changed, 24 insertions(+), 7 deletions(-)

diff --git a/modules/uncertainty/man/compute_sobol_indices.Rd b/modules/uncertainty/man/compute_sobol_indices.Rd
index 15fd8d350a8..b05e6c1f9a0 100644
--- a/modules/uncertainty/man/compute_sobol_indices.Rd
+++ b/modules/uncertainty/man/compute_sobol_indices.Rd
@@ -7,13 +7,13 @@
 compute_sobol_indices(outdir, sobol_obj, var = "GPP")
 }
 \arguments{
-\item{outdir}{PEcAn run output directory containing `ensemble.output.*.Rdata`
+\item{outdir}{PEcAn run output directory containing \code{ensemble.output.*.Rdata}
 files.}
 
 \item{sobol_obj}{object produced by
-`PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)`.}
+\code{PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)}.}
 
-\item{var}{Variable name to summarize (default `"GPP"`).}
+\item{var}{Variable name to summarize (default \code{"GPP"}).}
 }
 \value{
 A tibble of Sobol first-order and total-order indices with attached
@@ -21,15 +21,32 @@ A tibble of Sobol first-order and total-order indices with attached
 }
 \description{
 Loads standardized ensemble output for a Sobol run, computes first-order and
-total-order Sobol indices with `sensobol`, and saves both the Sobol design
+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
-parameter alone, while total-order indices summarize the full contribution of
-that parameter including its interactions with other parameters. PEcAn uses
-`sensobol` for these variance-based estimators; see Saltelli et al. (2008)
+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,

From e708bc075f47eb672dfc3a96a1815079697d04ae Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:18:16 -0400
Subject: [PATCH 22/25] update generate_joint_ensemble_design.Rd

---
 .../man/generate_joint_ensemble_design.Rd     | 20 ++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/modules/uncertainty/man/generate_joint_ensemble_design.Rd b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
index 145796d64ae..a1e3406bee7 100644
--- a/modules/uncertainty/man/generate_joint_ensemble_design.Rd
+++ b/modules/uncertainty/man/generate_joint_ensemble_design.Rd
@@ -2,10 +2,7 @@
 % 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)
 }
@@ -14,7 +11,8 @@ generate_joint_ensemble_design(settings, ensemble_size, sobol = FALSE)
 
 \item{ensemble_size}{Integer specifying the number of ensemble members.
 When \code{sobol = TRUE}, this is the Sobol base sample size \code{N}, not
-the expanded number of model runs.}
+the expanded number of model runs (which will be \code{N * (k + 2)} for
+\code{k} independent factors).}
 
 \item{sobol}{Logical, generate a variance-based Sobol design using
 \code{sensobol}.}
@@ -23,11 +21,19 @@ the expanded number of model runs.}
 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{\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{
+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.
+}

From 2d253ab916b2f79e830ff2d25736d4650a515192 Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 11:18:39 -0400
Subject: [PATCH 23/25] add trait_sample_bank_size.Rd

---
 .../uncertainty/man/trait_sample_bank_size.Rd | 21 +++++++++++++++++++
 1 file changed, 21 insertions(+)
 create mode 100644 modules/uncertainty/man/trait_sample_bank_size.Rd

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}

From fd4289c589b1229381aec9c955ee57522324c8eb Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 15:07:39 -0400
Subject: [PATCH 24/25] update NEWS.md

---
 modules/uncertainty/NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/modules/uncertainty/NEWS.md b/modules/uncertainty/NEWS.md
index 92cab17f25c..4614fc8e83a 100644
--- a/modules/uncertainty/NEWS.md
+++ b/modules/uncertainty/NEWS.md
@@ -1,7 +1,7 @@
 # PEcAn.uncertainty 1.9.0.9000
 
 * Switched Sobol global sensitivity workflows to `sensobol`, including explicit design metadata and saved Sobol index outputs.
-* Fixed parent-linked input design handling and validation so parameter-selected sample indices stay aligned and in bounds during workflow configuration.
+* 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

From ae29eba4633338113e04a91021c3c9618b6c400a Mon Sep 17 00:00:00 2001
From: divne7022 <www.divinehome12@gmail.com>
Date: Tue, 24 Mar 2026 15:08:02 -0400
Subject: [PATCH 25/25] update CHANGELOG.md

---
 CHANGELOG.md | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2f770868d06..83f2c3f61a9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,7 +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 `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).
@@ -19,9 +19,6 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Added datasets to `PEcAn.data.land` 
   *  `landiq_crop_mapping_codes` dataset mapping LandIQ crop classification codes to human-readable crop names.
   *  `bism_kc_by_crop` dataset containing BISm crop coefficient schedules and stage timing references for use in ET estimation, including columns that map to LandIQ class and subclass.
-
-### Fixed
-- Fixed joint ensemble design and workflow configuration handling to preserve parent-child input alignment and validate `input_design$param` indices before writing configs.
 - PEcAn.SIPNET gains support for SIPNET v2, whose features includes management events, nitrogen cycle tracking, explicit N2O and methane fluxes, runtime setting of feature flags, and changes to the parameter set (now 73 parameters). SIPNET v1 is still fully supported, but workarounds for bugs in the legacy `sipnet.unk` version have been removed.
 
 ### Fixed