diff --git a/.Rbuildignore b/.Rbuildignore
index 739ce55..2697148 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -6,3 +6,5 @@
 ^\.github$
 ^data-raw$
 ^README\.Rmd$
+^doc$
+^Meta$
diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml
index 865bcf6..6d574db 100644
--- a/.github/workflows/R-CMD-check.yml
+++ b/.github/workflows/R-CMD-check.yml
@@ -26,7 +26,7 @@ jobs:
           - {os: windows-latest,  r: 'release'}
           - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
           - {os: ubuntu-latest,   r: 'release'}
-          - {os: ubuntu-20.04,   r: '3.6'}
+          - {os: ubuntu-22.04,    r: '3.6'}
           
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
@@ -42,7 +42,7 @@ jobs:
 
       - uses: actions/setup-python@v4
         with:
-          python-version: '3.11' # Version range or exact version of a Python version to use, using SemVer's version range syntax
+          python-version: '3.14' # Version range or exact version of a Python version to use, using SemVer's version range syntax
           architecture: 'x64' # optional x64 or x86. Defaults to x64 if not specified
 
       - uses: r-lib/actions/setup-pandoc@v2
diff --git a/.gitignore b/.gitignore
index f6345fa..1acb996 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,6 @@
 .RData
 .Ruserdata
 ^rosettaPTF\.Rproj$
+README.html
+/doc/
+/Meta/
diff --git a/DESCRIPTION b/DESCRIPTION
index 811c6e4..5407b0f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,9 +1,9 @@
 Package: rosettaPTF
 Title: R Frontend for Rosetta Pedotransfer Functions
-Version: 0.1.5
+Version: 0.2.0
 Author: Soil and Plant Science Division Staff
 Maintainer: Andrew G. Brown <andrew.g.brown@usda.gov>
-Description: Access Python rosetta-soil pedotransfer functions in an R environment. Rosetta is a neural network-based model for predicting unsaturated soil hydraulic parameters from basic soil characterization data. The model predicts parameters for the van Genuchten unsaturated soil hydraulic properties model, using sand, silt, and clay, bulk density and water content. The codebase is now maintained by Dr. Todd Skaggs and other U.S. Department of Agriculture employees. This R package is intended to provide for use cases that involve many thousands of calls to the pedotransfer function. Less demanding use cases are encouraged to use the web interface or API endpoint. There are additional wrappers of the API endpoints provided by the soilDB R package `ROSETTA()` method.
+Description: Access the rosetta-soil Python pedotransfer functions from R. Rosetta is a neural network-based model for predicting unsaturated soil hydraulic parameters from basic soil characterization data (sand, silt, clay, bulk density, and water content). Predictions are made for the van Genuchten unsaturated hydraulic properties model, with uncertainty quantification via bootstrap ensemble. Designed for efficient batch processing of large datasets through vectorized computation and optional parallel processing.
 Config/reticulate:
   list(
     packages = list(
@@ -15,7 +15,7 @@ License: GPL (>= 2)
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Depends: R (>= 3.5)
 URL: https://github.com/ncss-tech/rosettaPTF, https://ncss-tech.github.io/rosettaPTF/
 BugReports: https://github.com/ncss-tech/rosettaPTF/issues
@@ -24,4 +24,6 @@ Imports:
     reticulate,
     terra
 Suggests: 
+    litedown,
     testthat
+VignetteBuilder: litedown
diff --git a/NAMESPACE b/NAMESPACE
index 4b14cd5..a85d8f9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(ann_predict,Rosetta)
 S3method(ann_predict,default)
 S3method(predict,Rosetta)
+S3method(predict,UnsaturatedK)
 S3method(py_to_r,rosetta.rosetta.SoilData)
 S3method(run_rosetta,RasterBrick)
 S3method(run_rosetta,RasterStack)
@@ -11,10 +13,12 @@ S3method(run_rosetta,default)
 S3method(run_rosetta,matrix)
 export(Rosetta)
 export(SoilDataFromArray)
+export(UnsaturatedK)
 export(ann_predict)
 export(find_python)
 export(get_rosetta_module)
 export(install_rosetta)
+export(rosesoil)
 export(rosetta_module_available)
 export(run_rosetta)
 importFrom(parallel,makeCluster)
@@ -32,6 +36,7 @@ importFrom(reticulate,r_to_py)
 importFrom(reticulate,use_condaenv)
 importFrom(reticulate,use_python)
 importFrom(stats,na.omit)
+importFrom(stats,predict)
 importFrom(terra,`nlyr<-`)
 importFrom(terra,rast)
 importFrom(terra,readStart)
diff --git a/NEWS.md b/NEWS.md
index e143a61..078d09a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,14 @@
+# rosettaPTF 0.2.0
+
+* Compatible with `rosetta-soil` Python package v0.3
+  * Updated `run_rosetta()` and `predict.Rosetta()` to handle the 7-parameter output (adding `K0` and `L`) introduced in `rosetta-soil` v0.3.
+  * Added `estimate_type` argument to `run_rosetta()` to support linear, logarithmic (default), and geometric parameter estimations. Improved documentation and added examples for the geometric scale.
+  * Added `UnsaturatedK()` R constructor and `predict.UnsaturatedK()` method for predicting `K0` and `L` from retention parameters (requires `rosetta-soil` >= 0.3).
+  * Added `rosesoil()` R wrapper for the new upstream `rosesoil()` function (requires `rosetta-soil` >= 0.3).
+  * Deprecated `SoilDataFromArray()` in favor of direct list input (supported in `rosetta-soil` >= 0.3.
+  * Deprecated `ann_predict()` as the underlying Python method has been removed in v0.3. It now redirects to `predict()`.
+* Added a new vignette: **"Performance Optimization and Raster Processing"** covering best practices for high-throughput workflows.
+
 # rosettaPTF 0.1.5
 
 * Fix check logic for whether input SpatRaster is in memory
diff --git a/R/AAAA.R b/R/AAAA.R
index 07fc060..82353cb 100644
--- a/R/AAAA.R
+++ b/R/AAAA.R
@@ -31,7 +31,6 @@ numpy_module <- NULL
   }
 
   !is.null(rosetta_module) && !is.null(numpy_module)
-
 }
 
 #' @importFrom reticulate configure_environment
diff --git a/R/Class-Rosetta.R b/R/Class-Rosetta.R
index 5cf54c3..83a5846 100644
--- a/R/Class-Rosetta.R
+++ b/R/Class-Rosetta.R
@@ -16,7 +16,7 @@
 #' @rdname Rosetta-class
 #' @export
 Rosetta <- function(rosetta_version = 3, model_code = 3) {
-  object <- rosetta_module$Rosetta(rosetta_version, model_code)
+  object <- rosetta_module$Rosetta(as.integer(rosetta_version), as.integer(model_code))
   structure(list(object = object), class = "Rosetta")
 }
 
@@ -24,14 +24,38 @@ Rosetta <- function(rosetta_version = 3, model_code = 3) {
 #' @param object _Rosetta_ object containing class instance (e.g. from `Rosetta()`)
 #' @param soildata A list containing vectors; with number of parameters matching the model type of `object`
 #' @param ... not used
+#' @return A list containing `mean` and `stdev` matrices (one row per sample).
+#'
+#' For `rosetta-soil` >= 0.3, the columns are: `theta_r`, `theta_s`, `alpha`, `npar`, `ksat`.
+#' Note that these parameters are in the scale produced by the underlying model (often log10 for alpha, npar, and ksat).
 #' @importFrom reticulate r_to_py import
 #' @method predict Rosetta
 #' @export
 #' @examples
 #' # predict(Rosetta(), list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
 predict.Rosetta <- function(object, soildata, ...) {
-  object$object$predict(numpy_module$array(reticulate::r_to_py(soildata),
-                                           dtype = "float"))
+  if (rosetta_pkg_version() >= package_version("0.3.0")) {
+    res <- object$object$predict(numpy_module$array(reticulate::r_to_py(soildata),
+                                             dtype = "float"))
+    retc_boot <- res[[1]]
+    ksat_boot <- res[[2]]
+
+    retc_mean <- numpy_module$mean(retc_boot, axis = 0L)
+    retc_std  <- numpy_module$std(retc_boot, axis = 0L)
+    ksat_mean <- numpy_module$mean(ksat_boot, axis = 0L)
+    ksat_std  <- numpy_module$std(ksat_boot, axis = 0L)
+
+    mean_val <- numpy_module$concatenate(list(retc_mean, ksat_mean), axis = 1L)
+    std_val  <- numpy_module$concatenate(list(retc_std, ksat_std), axis = 1L)
+
+    return(list(mean = mean_val, stdev = std_val))
+
+  } else {
+    res <- object$object$predict(numpy_module$array(reticulate::r_to_py(soildata),
+                                             dtype = "float"))
+    names(res) <- c("mean", "stdev")
+    return(res)
+  }
 }
 
 #' Extended _Rosetta_ Predictions, Parameter Distributions and Summary Statistics after Zhang & Schaap (2017)
@@ -47,16 +71,65 @@ ann_predict <- function(object, soildata, sum_data = TRUE)
 #' @rdname ann_predict
 #' @export
 ann_predict.default <- function(object, soildata, sum_data = TRUE) {
-  message("ann_predict() is defined for objects with class Rosetta; see `Rosetta()` to create a new instance")
+  if (rosetta_pkg_version() >= package_version("0.3.0")) {
+    .Deprecated("predict", msg = "ann_predict() is deprecated in rosetta-soil >= 0.3.0. Use predict() instead.")
+  } else {
+    message("ann_predict() is defined for objects with class Rosetta; see `Rosetta()` to create a new instance")
+  }
   ann_predict.Rosetta(object = object, soildata = soildata, sum_data = sum_data)
 }
 
 #' @rdname ann_predict
 #' @method ann_predict Rosetta
+#' @export
+#' @importFrom stats predict
 #' @examples
 #' # ann_predict(Rosetta(), list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
 ann_predict.Rosetta <- function(object, soildata, sum_data = TRUE) {
+  if (rosetta_pkg_version() >= package_version("0.3.0")) {
+    .Deprecated("predict", msg = "ann_predict() is deprecated in rosetta-soil >= 0.3.0. Use predict() instead.")
+    return(predict(object, soildata))
+  }
   object$object$ann_predict(numpy_module$array(reticulate::r_to_py(soildata),
                                                dtype = "float"),
                             sum_data = sum_data)
 }
+
+#' Make an UnsaturatedK object instance
+#'
+#' @description `UnsaturatedK`: Create an instance of the `UnsaturatedK` class from `rosetta-soil` >= 0.3. This class is used to predict `K0` and `L` from retention parameters.
+#'
+#' @return an instance of the `UnsaturatedK` class.
+#' @export
+UnsaturatedK <- function() {
+  if (rosetta_pkg_version() < package_version("0.3.0")) {
+    stop("UnsaturatedK requires rosetta-soil >= 0.3.0", call. = FALSE)
+  }
+
+  object <- rosetta_module$UnsaturatedK()
+  structure(list(object = object), class = "UnsaturatedK")
+}
+
+#' Predict K0 and L from retention parameters
+#'
+#' @param object _UnsaturatedK_ object
+#' @param retc_params A list or matrix of retention parameters (theta_r, theta_s, alpha, npar)
+#' @param ... not used
+#' @return a `data.frame` with `log10_K0_mean`, `lpar_mean`, `log10_K0_sd`, `lpar_sd`
+#' @method predict UnsaturatedK
+#' @export
+predict.UnsaturatedK <- function(object, retc_params, ...) {
+  res <- object$object$predict(numpy_module$array(reticulate::r_to_py(retc_params),
+                                           dtype = "float"))
+
+  k0l_mean <- numpy_module$mean(res, axis = 0L)
+  k0l_std  <- numpy_module$std(res, axis = 0L)
+
+  df <- data.frame(
+    log10_K0_mean = k0l_mean[, 1],
+    lpar_mean     = k0l_mean[, 2],
+    log10_K0_sd   = k0l_std[, 1],
+    lpar_sd       = k0l_std[, 2]
+  )
+  return(df)
+}
diff --git a/R/rosesoil.R b/R/rosesoil.R
new file mode 100644
index 0000000..5094e88
--- /dev/null
+++ b/R/rosesoil.R
@@ -0,0 +1,58 @@
+#' Run rosesoil() from rosetta-soil >= 0.3.0
+#'
+#' @param soildata A list of numeric vectors or a data.frame (3-6 columns: sand, silt, clay, optionally bulk density, th33, and th1500)
+#' @param rosetta_version integer, 1-3. Default: 3
+#' @param estimate_type _character_. One of `"arith"` (default), `"log"`, or `"geo"`. Only used if `rosetta-soil` >= 0.3.1. `"log"` returns parameters on a logarithmic (log10) scale for `alpha`, `npar`, `ksat`, and `k0`. `"geo"` returns the geometric mean of bootstrap estimates (exponent of the mean of log-transformed values). This is often preferred for parameters that vary by orders of magnitude, such as `alpha` and `ksat`.
+#' @param vars optional column name mapping (same as run_rosetta)
+#' @return a data.frame with all RosettaResult fields
+#' @export
+rosesoil <- function(soildata, rosetta_version = 3, estimate_type = "arith", vars = NULL) {
+  if (rosetta_pkg_version() < package_version("0.3.0")) {
+    stop("rosesoil() requires rosetta-soil >= 0.3.0. Please run install_rosetta(upgrade = TRUE).")
+  }
+
+  if (inherits(soildata, "data.frame")) {
+    if (!is.null(vars)) {
+      if (!all(vars %in% colnames(soildata))) {
+        stop("all custom parameter names in `vars` must be present in `soildata`",
+             call. = FALSE)
+      } else {
+        soildata <- soildata[, vars[seq_along(colnames(soildata))]]
+      }
+    }
+
+    nid <- nrow(soildata)
+    soildatatemplate <- data.frame(
+      sand = numeric(nid),
+      silt = numeric(nid),
+      clay = numeric(nid),
+      bulkdensity = numeric(nid),
+      th33 = numeric(nid),
+      th1500 = numeric(nid)
+    )
+    soildatatemplate[] <- NA_real_
+    soildatatemplate[, 1:ncol(soildata)] <- soildata
+    soildata_list <- unlist(apply(soildatatemplate, 1,
+                               function(x)
+                                 list(as.numeric(
+                                   stats::na.omit(as.numeric(x))
+                                 ))),
+                         recursive = FALSE)
+  } else {
+    soildata_list <- soildata
+  }
+
+  res_obj <- rosetta_module$rosesoil(as.integer(rosetta_version),
+                                    soildata_list,
+                                    estimate_type = estimate_type)
+
+  res_dicts <- res_obj$asdicts()
+
+  # handle NULL values in dicts (convert to NA)
+  res_df <- do.call(rbind, lapply(res_dicts, function(d) {
+    d[sapply(d, is.null)] <- NA_real_
+    as.data.frame(d)
+  }))
+
+  return(res_df)
+}
diff --git a/R/rosetta_utils.R b/R/rosetta_utils.R
index 3f5900d..9792e20 100644
--- a/R/rosetta_utils.R
+++ b/R/rosetta_utils.R
@@ -2,12 +2,16 @@
 
 #' Convert list of numeric vectors to _SoilData_ Python object
 #'
-#' @description `SoilDataFromArray`: convert a list of numeric vectors containing soil properties to a `rosetta.rosetta.SoilData` class
+#' @description `SoilDataFromArray`: convert a list of numeric vectors containing soil properties to a `rosetta.rosetta.SoilData` class. In `rosetta-soil` >= 0.3, direct list input is preferred.
 #'
 #' @param x a list of numeric vectors
 #' @return an object reference to a Rosetta _SoilData_ Python object constructed from `x`
 #' @export
 SoilDataFromArray <- function(x) {
+  if (rosetta_pkg_version() >= package_version("0.3.0")) {
+    .Deprecated(msg = "Direct list input is now supported by rosetta-soil >= 0.3.0. SoilDataFromArray is deprecated.")
+    return(x)
+  }
   rosetta_module$SoilData$from_array(x)
 }
 
@@ -20,6 +24,20 @@ py_to_r.rosetta.rosetta.SoilData <- function(x) {
   x
 }
 
+#' Get rosetta-soil Python package version
+#' @return `package_version` object
+#' @keywords internal
+rosetta_pkg_version <- function() {
+  if (rosetta_module_available()) {
+    v <- try(rosetta_module$`__version__`, silent = TRUE)
+    if (inherits(v, "try-error") || is.null(v)) {
+      return(package_version("0.1.0"))
+    }
+    return(package_version(v))
+  }
+  package_version("0.0.0")
+}
+
 #' Check if Rosetta module is available for import from local Python environment
 #' @return _logical_
 #' @export
diff --git a/R/run_rosetta.R b/R/run_rosetta.R
index 9ae541d..3c3e87b 100644
--- a/R/run_rosetta.R
+++ b/R/run_rosetta.R
@@ -3,40 +3,66 @@
 #' @param soildata A list of numeric vectors each containing 3 to 6 values: `"sand"`, `"silt"`, `"clay"`, `"bulkdensity"`, `"th33"`, `"th1500"`, a _data.frame_ or _matrix_ with 3 to 6 columns OR a `Raster*`/`SpatRaster` object with 3 to 6 layers. Sand, silt, and clay must sum to a total of 100%.
 #' @param vars _character_. Optional: names and order of custom column names if `soildata` is a _data.frame_, _RasterStack_, _RasterBrick_ or _SpatRaster_. Default `NULL` assumes input column order follows `sand`, `silt`, `clay`, `bulkdensity`, `th33`, `th1500` and does not check names.
 #' @param rosetta_version Default: 3
+#' @param estimate_type _character_. One of `"log"` (default), `"arith"`, or `"geo"`. Only used if `rosetta-soil` >= 0.3.1. Default `"log"` preserves logarithmic (log10) scale for `alpha`, `npar`, and `Ksat`. `"geo"` returns the geometric mean of bootstrap estimates (exponent of the mean of log-transformed values). This is often preferred for parameters that vary by orders of magnitude, such as `alpha` and `Ksat`.
 #' @param ... additional arguments not used
 #'
-#' @return A _data.frame_ containing `mean` and `stdev` for following five columns (parameters for van Genuchten-Mualem equation)
+#' @return A _data.frame_ containing `mean` and `stdev` for the following columns (parameters for van Genuchten-Mualem equation)
 #' -	`"theta_r"`, residual water content
 #' -	`"theta_s"`, saturated water content
-#' -	`"log10(alpha)"`, 'alpha' shape parameter, log10(1/cm)
-#' -	`"log10(npar)"`, 'n' shape parameter
-#' -	`"log10(Ksat)"`, saturated hydraulic conductivity, log10(cm/day)
+#' -	`"alpha"`, 'alpha' shape parameter (1/cm). Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+#' -	`"npar"`, 'n' shape parameter. Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+#' -	`"Ksat"`, saturated hydraulic conductivity (cm/day). Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+#' -  `"K0"`, unsaturated hydraulic conductivity (cm/day). Only if `rosetta-soil` >= 0.3.1. Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+#' -  `"lpar"`, unsaturated hydraulic conductivity exponent. Only if `rosetta-soil` >= 0.3.1.
 #'
 #' If the sum of sand, silt, and clay is not 100%, the parameter value estimates will be `NaN`.
 #'
+#' @details
+#' ## Performance Note
+#'
+#' Use `cores > 1` with `SpatRaster` or `Raster*` inputs to parallelize processing of cells across multiple cores.
+#'
 #' @aliases run_rosetta
 #' @rdname run_rosetta
 #' @export
 run_rosetta.default <- function(soildata,
                                 vars = NULL,
-                                rosetta_version = 3, ...) {
+                                rosetta_version = 3,
+                                estimate_type = "log", ...) {
 
   if (is.numeric(soildata)) {
     soildata <- as.data.frame(t(soildata))
-    run_rosetta.data.frame(soildata = soildata, vars = vars, rosetta_version = rosetta_version)
+    return(run_rosetta.data.frame(soildata = soildata, vars = vars, rosetta_version = rosetta_version, estimate_type = estimate_type))
   }
 
   # identify records with enough data
   good.idx <- which(sapply(soildata, length) >= 3)
 
-  # run rosetta
-  res <- rosetta_module$rosetta(rosetta_version, SoilDataFromArray(soildata[good.idx]))
+  if (rosetta_pkg_version() >= package_version("0.3.1")) {
+    res <- rosetta_module$rosetta(as.integer(rosetta_version),
+                                 soildata[good.idx],
+                                 estimate_type = estimate_type)
+  } else {
+    res <- rosetta_module$rosetta(as.integer(rosetta_version), SoilDataFromArray(soildata[good.idx]))
+  }
 
   if (length(res) == 3) {
     names(res) <- c("mean","stdev","model_codes")
-    param_names <- c("theta_r", "theta_s", "log10_alpha", "log10_npar", "log10_Ksat")
+
+    nc <- ncol(res[[1]])
+    if (nc == 7) {
+       if (estimate_type == "log") {
+         param_names <- c("theta_r", "theta_s", "log10_alpha", "log10_npar", "log10_Ksat", "log10_K0", "lpar")
+       } else {
+         param_names <- c("theta_r", "theta_s", "alpha", "npar", "ksat", "k0", "lpar")
+       }
+    } else {
+      param_names <- c("theta_r", "theta_s", "log10_alpha", "log10_npar", "log10_Ksat")
+    }
+
     res[[1]] <- as.data.frame(res[[1]])
     colnames(res[[1]]) <- paste0(param_names, "_mean")
+    res[[2]] <- as.data.frame(res[[2]])
     colnames(res[[2]]) <- paste0(param_names, "_sd")
 
     res <- data.frame(model_code = res[[3]], cbind(res[[1]], res[[2]]))
@@ -56,6 +82,7 @@ run_rosetta.default <- function(soildata,
 run_rosetta <- function(soildata,
                         vars = NULL,
                         rosetta_version = 3,
+                        estimate_type = "log",
                         cores = 1,
                         core_thresh = NULL,
                         file = NULL,
@@ -69,6 +96,7 @@ run_rosetta <- function(soildata,
 run_rosetta.data.frame <- function(soildata,
                                    vars = NULL,
                                    rosetta_version = 3,
+                                   estimate_type = "log",
                                    ...) {
 
   # soildata <- as.data.frame(soildata)
@@ -101,13 +129,15 @@ run_rosetta.data.frame <- function(soildata,
   )
   soildatatemplate[] <- NA_real_
   soildatatemplate[, 1:ncol(soildata)] <- soildata
-  soildatain <- unlist(apply(soildatatemplate, 1,
-                             function(x)
-                               list(as.numeric(
-                                 na.omit(as.numeric(x))
-                               ))),
-                       recursive = FALSE)
-  run_rosetta.default(soildatain, vars = vars, rosetta_version = rosetta_version)
+
+  # Get number of non-NA columns per row
+  n_cols <- rowSums(!is.na(soildatatemplate))
+  m <- as.matrix(soildatatemplate)
+  soildatain <- lapply(seq_len(nid), function(i) {
+    m[i, 1:n_cols[i]]
+  })
+
+  run_rosetta.default(soildatain, vars = vars, rosetta_version = rosetta_version, estimate_type = estimate_type)
 }
 
 #' @export
@@ -115,8 +145,9 @@ run_rosetta.data.frame <- function(soildata,
 run_rosetta.matrix <- function(soildata,
                                vars = NULL,
                                rosetta_version = 3,
+                               estimate_type = "log",
                                ...) {
-  run_rosetta(as.data.frame(soildata), vars = vars, rosetta_version = 3)
+  run_rosetta(as.data.frame(soildata), vars = vars, rosetta_version = rosetta_version, estimate_type = estimate_type)
 }
 
 #' @export
@@ -125,25 +156,17 @@ run_rosetta.matrix <- function(soildata,
 run_rosetta.RasterStack <- function(soildata,
                                     vars = NULL,
                                     rosetta_version = 3,
+                                    estimate_type = "log",
                                     cores = 1,
                                     core_thresh = 20000L,
                                     file = paste0(tempfile(), ".tif"),
                                     nrows = nrow(soildata) / (terra::ncell(soildata) / core_thresh),
                                     overwrite = TRUE) {
-  ## for in memory only, can just convert to data.frame and use that method
-  # res <- run_rosetta(raster::as.data.frame(soildata),
-  #                    vars = vars,
-  #                    rosetta_version = rosetta_version)
-  # resstackout <- soildata
-  # for(i in 1:ncol(res)) {
-  #   resstackout[[i]] <- res[[i]]
-  # }
-  # names(resstackout) <- colnames(res)
-  # resstackout
   run_rosetta(
     terra::rast(soildata),
     vars = vars,
     rosetta_version = rosetta_version,
+    estimate_type = estimate_type,
     cores = cores,
     file = file,
     nrows = nrows,
@@ -157,6 +180,7 @@ run_rosetta.RasterStack <- function(soildata,
 run_rosetta.RasterBrick <- function(soildata,
                                     vars = NULL,
                                     rosetta_version = 3,
+                                    estimate_type = "log",
                                     cores = 1,
                                     core_thresh = 20000L,
                                     file = paste0(tempfile(), ".tif"),
@@ -165,6 +189,7 @@ run_rosetta.RasterBrick <- function(soildata,
   run_rosetta(terra::rast(soildata),
               vars = vars,
               rosetta_version = rosetta_version,
+              estimate_type = estimate_type,
               cores = cores,
               file = file,
               nrows = nrows,
@@ -183,6 +208,7 @@ run_rosetta.RasterBrick <- function(soildata,
 run_rosetta.SpatRaster <- function(soildata,
                                    vars = NULL,
                                    rosetta_version = 3,
+                                   estimate_type = "log",
                                    cores = 1,
                                    core_thresh = 20000L,
                                    file = paste0(tempfile(), ".tif"),
@@ -200,9 +226,11 @@ run_rosetta.SpatRaster <- function(soildata,
 
   # create template brick
   out <- terra::rast(soildata)
-  cnm <- c("id", "model_code", "theta_r_mean", "theta_s_mean", "log10_alpha_mean",
-           "log10_npar_mean", "log10_Ksat_mean", "theta_r_sd", "theta_s_sd",
-           "log10_alpha_sd", "log10_npar_sd", "log10_Ksat_sd")
+
+  # determine output columns by running a small sample
+  sample_res <- run_rosetta.default(list(c(33, 33, 34)), rosetta_version = rosetta_version, estimate_type = estimate_type)
+  cnm <- colnames(sample_res)
+
   terra::nlyr(out) <- length(cnm)
   names(out) <- cnm
   out_info <- terra::writeStart(out, filename = file, overwrite = overwrite)
@@ -220,20 +248,17 @@ run_rosetta.SpatRaster <- function(soildata,
     cls <- parallel::makeCluster(cores)
     on.exit(parallel::stopCluster(cls))
 
-    # TODO: can blocks be parallelized?
     for (i in seq_along(start_row)) {
       if (n_row[i] > 0) {
         blockdata <- terra::readValues(soildata, row = start_row[i], nrows = n_row[i], dataframe = TRUE)
-        ids <- 1:nrow(blockdata)
-        # soilDB makeChunks logic; what is tradeoff between chunk size and number of requests?
-        # run_rosetta is a "costly" function and not particularly fast, so in theory parallel would help
 
         # parallel within-block processing
-        n <- floor(length(ids) / core_thresh / cores) + 1
-        X <- split(blockdata, rep(seq(from = 1, to = n)))[1:length(ids)]
+        n <- max(cores, ceiling(nrow(blockdata) / core_thresh))
+        X <- split(blockdata, rep(seq_len(n), length.out = nrow(blockdata)))
         r <- do.call('rbind', parallel::clusterApply(cls, X, function(x) rosettaPTF::run_rosetta(x,
                                                                                                  vars = vars,
-                                                                                                 rosetta_version = rosetta_version)))
+                                                                                                 rosetta_version = rosetta_version,
+                                                                                                 estimate_type = estimate_type)))
 
         terra::writeValues(out, as.matrix(r), start_row[i], nrows = n_row[i])
       }
@@ -243,7 +268,8 @@ run_rosetta.SpatRaster <- function(soildata,
       if (n_row[i] > 0) {
         foo <- rosettaPTF::run_rosetta(terra::readValues(soildata, row = start_row[i], nrows = n_row[i], dataframe = TRUE),
                                        vars = vars,
-                                       rosetta_version = rosetta_version)
+                                       rosetta_version = rosetta_version,
+                                       estimate_type = estimate_type)
         terra::writeValues(out, as.matrix(foo), start_row[i], nrows = n_row[i])
       }
     }
diff --git a/README.Rmd b/README.Rmd
index e0888f9..c1f931c 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -1,11 +1,23 @@
 ---
 output: github_document
+knit: (function(input, ...) {
+    litedown::fuse(input, 'README.md')
+    x = readLines('README.md')
+    if (length(x) > 0 && x[1] == '---') {
+      i = grep('^---$', x)
+      if (length(i) >= 2) x = x[-(1:i[2])]
+    }
+    writeLines(x, 'README.md')
+  })
 ---
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-```{r, include = FALSE}
-knitr::opts_chunk$set(
+```{r setup, include = FALSE}
+library(rosettaPTF)
+library(terra)
+
+litedown::reactor(
   collapse = TRUE,
   comment = "#>",
   fig.path = "man/figures/README-",
@@ -46,13 +58,15 @@ library(rosettaPTF)
 
 The [rosetta-soil](https://github.com/usda-ars-ussl/rosetta-soil) module is a Python package maintained by Dr. Todd Skaggs (USDA-ARS) and other U.S. Department of Agriculture employees. 
 
-The Rosetta pedotransfer function predicts five parameters for the van Genuchten model of unsaturated soil hydraulic properties
+The Rosetta pedotransfer function predicts seven parameters (five in versions < 0.2.0) for the van Genuchten model of unsaturated soil hydraulic properties:
 
- - `theta_r` : residual volumetric water content
- - `theta_s` : saturated volumetric water content
- - `log10(alpha)` : retention shape parameter `[log10(1/cm)]`
- - `log10(n)` : retention shape parameter (also referred to as `npar`)
- - `log10(ksat)` : saturated hydraulic conductivity `[log10(cm/d)]`
+* `theta_r` : residual volumetric water content
+* `theta_s` : saturated volumetric water content
+* `alpha` : retention shape parameter `[1/cm]`. Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+* `npar` : retention shape parameter (also referred to as `n`). Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+* `ksat` : saturated hydraulic conductivity `[cm/d]`. Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+* `K0` : unsaturated hydraulic conductivity matching point `[cm/d]`. Logarithmic (log10) scale if `estimate_type="log"` (default); Geometric mean if `estimate_type="geo"`.
+* `lpar` : unsaturated hydraulic conductivity exponent.
 
 For each set of input data a mean and standard deviation of each parameter is given.
 
@@ -62,13 +76,10 @@ Less demanding use cases are encouraged to use the web interface or API endpoint
 
 The [Rosetta](http://ncss-tech.github.io/AQP/soilDB/ROSETTA-API.html) model relies on a minimum of 3 soil properties, with increasing (expected) accuracy as additional properties are included:
 
-  *  Required, `sand`, `silt`, `clay`: USDA soil texture separates (percentages) that sum to 100%
-
-  *  Optional, `bulk density (any moisture basis)`: mass per volume after accounting for >2mm fragments, units of grams/cm3
-
-  *  Optional, `volumetric water content at 33 kPa`: roughly “field capacity” for most soils, units of cm3/cm3
-
-  *  Optional, `volumetric water content at 1500 kPa`: roughly “permanent wilting point” for most plants, units of cm3/cm3
+  * Required, `sand`, `silt`, `clay`: USDA soil texture separates (percentages) that sum to 100%
+  * Optional, `bulk density (any moisture basis)`: mass per volume after accounting for >2mm fragments, units of grams/cm3
+  * Optional, `volumetric water content at 33 kPa`: roughly “field capacity” for most soils, units of cm3/cm3
+  * Optional, `volumetric water content at 1500 kPa`: roughly “permanent wilting point” for most plants, units of cm3/cm3
 
 The default order of inputs is: `sand`, `silt`, `clay`, `bulk density (any basis)`, `water content (field capacity; 33 kPa)`, `water content (permanent wilting point; 1500 kPa)` of which the first three are required. 
 
@@ -124,22 +135,48 @@ Alternately, to install the module manually with `pip` you can run the following
 python -m pip install rosetta-soil
 ```
 
+## High-Throughput Processing
+
+`{rosettaPTF}` supports efficient batch processing of large soil datasets through vectorized computation in the underlying `rosetta-soil` backend.
+
+For large datasets:
+
+* Use `cores > 1` with `run_rosetta()` and `SpatRaster` or `Raster*` inputs to parallelize the calls.
+
 ## `run_rosetta()`
 
-Batch runs of Rosetta models can be done using using `list`, `data.frame`, `matrix`, `RasterStack`, `RasterBrick` and `SpatRaster` objects as input. 
+Batch runs of Rosetta models can be done using `list`, `data.frame`, `matrix`, `RasterStack`, `RasterBrick` and `SpatRaster` objects as input.
+
+Plain R lists are the preferred input format. The helper `SoilDataFromArray()` is deprecated.
 
 ### `list()` Input Example
 
 ```{r}
+# Plain R lists are passed directly to Python
 run_rosetta(list(c(30, 30, 40, 1.5), c(55, 25, 20),  c(55, 25, 20, 1.1)),
             rosetta_version = 3)
 ```
 
 Output `model_code` reflects the number of parameters in the input. 
 
-### `data.frame()` Input Example
+### Parameter Estimation Scales
+
+By default, `{rosettaPTF}` uses `estimate_type = "log"` to maintain backward compatibility with previous versions, returning `alpha`, `npar`, and `Ksat` on a logarithmic (log10) scale. You can now request estimates on a linear scale directly:
+
+```{r}
+run_rosetta(list(c(30, 30, 40, 1.5)), estimate_type = "arith")
+```
+
+Note that the output column names will change to reflect the linear scale (e.g., `ksat_mean` instead of `log10_Ksat_mean`).
 
-The `data.frame` interface allows for using using custom column names and order. If the `vars` argument is not specified it is assumed that the columns are in the order specified in the `run_rosetta()` manual page.
+Additionally, `estimate_type = "geo"` can be used to return the **geometric mean** of the bootstrap estimates. This is often preferred for parameters like $K_{sat}$ and $\alpha$ which can span several orders of magnitude, as the geometric mean is less sensitive to extreme outliers in the bootstrap ensemble than the arithmetic mean. Mathematically, the geometric mean is equivalent to the exponent of the mean of the log-transformed values.
+
+```{r}
+run_rosetta(list(c(30, 30, 40, 1.5)), estimate_type = "geo")
+```
+
+The `data.frame` interface allows for using using custom column names and order.
+ If the `vars` argument is not specified it is assumed that the columns are in the order specified in the `run_rosetta()` manual page.
 
 ```{r}
 run_rosetta(data.frame(
@@ -231,16 +268,35 @@ predict(my_rosetta, list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
 
 ### Extended _Rosetta_ Predictions, Parameter Distributions and Summary Statistics after Zhang & Schaap (2017) with `ann_predict()`
 
+`ann_predict()` is deprecated and redirects to `predict()`, as the underlying bootstrap data is now returned by `predict()` and summarized by R.
+
 ```{r}
 ann_predict(my_rosetta, list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
 ```
 
-## Selected References
+## New Features in `rosetta-soil` 0.3.0
+
+### `rosesoil()`
+
+`rosesoil()` is a new R wrapper for the upstream `rosesoil()` function, which returns a structured result including all model metadata.
 
-Three versions of the ROSETTA model are available, selected using `rosetta_version` argument.
+```{r}
+rosesoil(list(c(33, 33, 34, 1.5)))
+```
+
+### `UnsaturatedK()`
+
+`UnsaturatedK` provides a way to predict unsaturated hydraulic conductivity parameters `K0` and `lpar` from retention parameters.
 
-  - `rosetta_version` 1 - Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3-4): 163-176. doi: 10.1016/S0022-1694(01)00466-8.
+```{r}
+uk <- UnsaturatedK()
+predict(uk, list(c(0.12, 0.42, 0.008, 1.29)))
+```
+
+## Selected References
 
-  - `rosetta_version` 2 - Schaap, M.G., A. Nemes, and M.T. van Genuchten. 2004. Comparison of Models for Indirect Estimation of Water Retention and Available Water in Surface Soils. Vadose Zone Journal 3(4): 1455-1463. doi: 10.2136/vzj2004.1455.
+Three versions of the ROSETTA model are available, selected using `rosetta_version` argument:
 
-  - `rosetta_version` 3 - Zhang, Y., and M.G. Schaap. 2017. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3). Journal of Hydrology 547: 39-53. doi: 10.1016/j.jhydrol.2017.01.004.
+  - `rosetta_version` 1: Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3-4): 163-176. doi: 10.1016/S0022-1694(01)00466-8.
+  - `rosetta_version` 2: Schaap, M.G., A. Nemes, and M.T. van Genuchten. 2004. Comparison of Models for Indirect Estimation of Water Retention and Available Water in Surface Soils. Vadose Zone Journal 3(4): 1455-1463. doi: 10.2136/vzj2004.1455.
+  - `rosetta_version` 3: Zhang, Y., and M.G. Schaap. 2017. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3). Journal of Hydrology 547: 39-53. doi: 10.1016/j.jhydrol.2017.01.004. Version 3 includes predictions for unsaturated conductivity parameters `K0` and `lpar`.
diff --git a/README.md b/README.md
index dac471f..f4a58e0 100644
--- a/README.md
+++ b/README.md
@@ -46,8 +46,6 @@ will be notified.
 
 ``` r
 library(rosettaPTF)
-#> Error in system2(command = python, args = shQuote(script), stdout = TRUE,  : 
-#>   'CreateProcess' failed to run 'C:\Users\ANDREW~1.BRO\ONEDRI~1\DOCUME~1\VIRTUA~1\R-RETI~1\Scripts\python.exe "C:/Users/Andrew.G.Brown/AppData/Local/R/win-library/4.3/reticulate/config/config.py"' "C:/Users/Andrew.G.Brown/OneDrive - USDA/Documents/.virtualenvs/r-reticulate/Scripts/python.exe"
 ```
 
 ### `rosetta-soil` Python module
@@ -56,14 +54,25 @@ The [rosetta-soil](https://github.com/usda-ars-ussl/rosetta-soil) module
 is a Python package maintained by Dr. Todd Skaggs (USDA-ARS) and other
 U.S. Department of Agriculture employees.
 
-The Rosetta pedotransfer function predicts five parameters for the van
-Genuchten model of unsaturated soil hydraulic properties
+The Rosetta pedotransfer function predicts seven parameters (five in
+versions \< 0.2.0) for the van Genuchten model of unsaturated soil
+hydraulic properties:
 
 - `theta_r` : residual volumetric water content
 - `theta_s` : saturated volumetric water content
-- `log10(alpha)` : retention shape parameter `[log10(1/cm)]`
-- `log10(n)` : retention shape parameter (also referred to as `npar`)
-- `log10(ksat)` : saturated hydraulic conductivity `[log10(cm/d)]`
+- `alpha` : retention shape parameter `[1/cm]`. Logarithmic (log10)
+  scale if `estimate_type="log"` (default); Geometric mean if
+  `estimate_type="geo"`.
+- `npar` : retention shape parameter (also referred to as `n`).
+  Logarithmic (log10) scale if `estimate_type="log"` (default);
+  Geometric mean if `estimate_type="geo"`.
+- `ksat` : saturated hydraulic conductivity `[cm/d]`. Logarithmic
+  (log10) scale if `estimate_type="log"` (default); Geometric mean if
+  `estimate_type="geo"`.
+- `K0` : unsaturated hydraulic conductivity matching point `[cm/d]`.
+  Logarithmic (log10) scale if `estimate_type="log"` (default);
+  Geometric mean if `estimate_type="geo"`.
+- `lpar` : unsaturated hydraulic conductivity exponent.
 
 For each set of input data a mean and standard deviation of each
 parameter is given.
@@ -82,13 +91,10 @@ model relies on a minimum of 3 soil properties, with increasing
 
 - Required, `sand`, `silt`, `clay`: USDA soil texture separates
   (percentages) that sum to 100%
-
 - Optional, `bulk density (any moisture basis)`: mass per volume after
   accounting for \>2mm fragments, units of grams/cm3
-
 - Optional, `volumetric water content at 33 kPa`: roughly “field
   capacity” for most soils, units of cm3/cm3
-
 - Optional, `volumetric water content at 1500 kPa`: roughly “permanent
   wilting point” for most plants, units of cm3/cm3
 
@@ -123,9 +129,10 @@ reticulate::virtualenv_create("r-reticulate")
 
 ``` r
 rosettaPTF::find_python()
-#> [1] "C:/Program Files/Python312/python.exe"
 ```
 
+    ## [1] "/home/andrew/.virtualenvs/r-reticulate/bin/python"
+
 `find_python()` provides heuristics for setting up {reticulate} to use
 Python in commonly installed locations.
 
@@ -160,10 +167,14 @@ module you should restart your R session.
 
 ``` r
 rosettaPTF::install_rosetta()
-#> Using virtual environment "~/.virtualenvs/r-reticulate" ...
-#> [1] TRUE
 ```
 
+    ## Using virtual environment '/home/andrew/.virtualenvs/r-reticulate' ...
+
+    ## + /home/andrew/.virtualenvs/r-reticulate/bin/python -m pip install --upgrade --no-user NA --upgrade 'rosetta-soil==0.3.1'
+
+    ## [1] TRUE
+
 Alternately, to install the module manually with `pip` you can run the
 following command. This assumes a Python 3 binary called `python` can be
 found on your path.
@@ -172,34 +183,88 @@ found on your path.
 python -m pip install rosetta-soil
 ```
 
+## High-Throughput Processing
+
+`{rosettaPTF}` supports efficient batch processing of large soil
+datasets through vectorized computation in the underlying `rosetta-soil`
+backend.
+
+For large datasets:
+
+- Use `cores > 1` with `run_rosetta()` and `SpatRaster` or `Raster*`
+  inputs to parallelize the calls.
+
 ## `run_rosetta()`
 
-Batch runs of Rosetta models can be done using using `list`,
-`data.frame`, `matrix`, `RasterStack`, `RasterBrick` and `SpatRaster`
-objects as input.
+Batch runs of Rosetta models can be done using `list`, `data.frame`,
+`matrix`, `RasterStack`, `RasterBrick` and `SpatRaster` objects as
+input.
+
+Plain R lists are the preferred input format. The helper
+`SoilDataFromArray()` is deprecated.
 
 ### `list()` Input Example
 
 ``` r
+# Plain R lists are passed directly to Python
 run_rosetta(list(c(30, 30, 40, 1.5), c(55, 25, 20),  c(55, 25, 20, 1.1)),
             rosetta_version = 3)
-#>   id model_code theta_r_mean theta_s_mean log10_alpha_mean log10_npar_mean
-#> 1  1          3   0.11535773    0.4179120        -2.067139       0.1120102
-#> 2  2          2   0.08613275    0.3888528        -1.898150       0.1347136
-#> 3  3          3   0.09130753    0.4850320        -2.022388       0.1510716
-#>   log10_Ksat_mean  theta_r_sd  theta_s_sd log10_alpha_sd log10_npar_sd
-#> 1       0.8325407 0.013350113 0.009377977     0.08251142    0.01323413
-#> 2       1.1858005 0.006014445 0.006273536     0.07481303    0.01160419
-#> 3       1.9060148 0.012771407 0.013062171     0.10020312    0.01763982
-#>   log10_Ksat_sd
-#> 1    0.09245277
-#> 2    0.08428578
-#> 3    0.14163567
 ```
 
+    ##   id model_code theta_r_mean theta_s_mean log10_alpha_mean log10_npar_mean
+    ## 1  1          3   0.11535773    0.4179120        -2.067139       0.1120102
+    ## 2  2          2   0.08613275    0.3888528        -1.898150       0.1347136
+    ## 3  3          3   0.09130753    0.4850320        -2.022388       0.1510716
+    ##   log10_Ksat_mean log10_K0_mean  lpar_mean  theta_r_sd  theta_s_sd
+    ## 1       0.8325407    0.02327444 -1.0349735 0.013356794 0.009382669
+    ## 2       1.1858005    0.41950414 -0.9533283 0.006017454 0.006276675
+    ## 3       1.9060148    0.37302013 -0.3253938 0.012777798 0.013068707
+    ##   log10_alpha_sd log10_npar_sd log10_Ksat_sd log10_K0_sd  lpar_sd
+    ## 1     0.08255271    0.01324075    0.09249903   0.2317968 1.588173
+    ## 2     0.07485046    0.01161000    0.08432796   0.2275761 1.099102
+    ## 3     0.10025326    0.01764865    0.14170654   0.2540277 1.244062
+
 Output `model_code` reflects the number of parameters in the input.
 
-### `data.frame()` Input Example
+### Parameter Estimation Scales
+
+By default, `{rosettaPTF}` uses `estimate_type = "log"` to maintain
+backward compatibility with previous versions, returning `alpha`,
+`npar`, and `Ksat` on a logarithmic (log10) scale. You can now request
+estimates on a linear scale directly:
+
+``` r
+run_rosetta(list(c(30, 30, 40, 1.5)), estimate_type = "arith")
+```
+
+    ##   id model_code theta_r_mean theta_s_mean  alpha_mean npar_mean ksat_mean
+    ## 1  1          3    0.1153577     0.417912 0.008722012  1.294826  6.954245
+    ##    k0_mean lpar_mean theta_r_sd  theta_s_sd    alpha_sd    npar_sd  ksat_sd
+    ## 1 1.250061 -1.048504 0.01335679 0.009382669 0.001649266 0.03939482 1.470258
+    ##       k0_sd  lpar_sd
+    ## 1 0.7856297 1.576854
+
+Note that the output column names will change to reflect the linear
+scale (e.g., `ksat_mean` instead of `log10_Ksat_mean`).
+
+Additionally, `estimate_type = "geo"` can be used to return the
+**geometric mean** of the bootstrap estimates. This is often preferred
+for parameters like $K_{sat}$ and $\alpha$ which can span several orders
+of magnitude, as the geometric mean is less sensitive to extreme
+outliers in the bootstrap ensemble than the arithmetic mean.
+Mathematically, the geometric mean is equivalent to the exponent of the
+mean of the log-transformed values.
+
+``` r
+run_rosetta(list(c(30, 30, 40, 1.5)), estimate_type = "geo")
+```
+
+    ##   id model_code theta_r_mean theta_s_mean  alpha_mean npar_mean ksat_mean
+    ## 1  1          3    0.1153577     0.417912 0.008567645  1.294226  6.800498
+    ##    k0_mean lpar_mean theta_r_sd  theta_s_sd    alpha_sd    npar_sd  ksat_sd
+    ## 1 1.055053 -1.034973 0.01335679 0.009382669 0.001649266 0.03939482 1.470258
+    ##       k0_sd  lpar_sd
+    ## 1 0.7666369 1.588173
 
 The `data.frame` interface allows for using using custom column names
 and order. If the `vars` argument is not specified it is assumed that
@@ -213,17 +278,18 @@ run_rosetta(data.frame(
   a = 20,
   c = 20
 ), vars = letters[1:4])
-#>   id model_code theta_r_mean theta_s_mean log10_alpha_mean log10_npar_mean
-#> 1  1          2   0.08994502    0.4301366        -2.426236       0.1756873
-#> 2  2          3   0.08495731    0.3887858        -2.318826       0.1598879
-#>   log10_Ksat_mean  theta_r_sd  theta_s_sd log10_alpha_sd log10_npar_sd
-#> 1       1.1927311 0.006707593 0.008785824     0.07413139    0.01323068
-#> 2       0.9961317 0.010184683 0.008100061     0.07976954    0.01753829
-#>   log10_Ksat_sd
-#> 1    0.08709446
-#> 2    0.07771481
 ```
 
+    ##   id model_code theta_r_mean theta_s_mean log10_alpha_mean log10_npar_mean
+    ## 1  1          2   0.08994502    0.4301366        -2.426236       0.1756873
+    ## 2  2          3   0.08495731    0.3887858        -2.318826       0.1598879
+    ##   log10_Ksat_mean log10_K0_mean  lpar_mean  theta_r_sd  theta_s_sd
+    ## 1       1.1927311   -0.10923995  0.1813931 0.006710949 0.008790221
+    ## 2       0.9961317   -0.03714337 -0.1213272 0.010189780 0.008104114
+    ##   log10_alpha_sd log10_npar_sd log10_Ksat_sd log10_K0_sd  lpar_sd
+    ## 1     0.07416849    0.01323730    0.08713804   0.2194768 1.622813
+    ## 2     0.07980946    0.01754707    0.07775370   0.2056933 1.363478
+
 ### Soil Data Access / SSURGO Mapunit Aggregate Input Example
 
 This example pulls mapunit/component data from Soil Data Access (SDA).
@@ -236,37 +302,36 @@ results (1:1 with `mukey`).
 ``` r
 library(soilDB)
 library(terra)
-#> Warning: package 'terra' was built under R version 4.3.3
-#> terra 1.7.78
-```
-
-``` r
 library(rosettaPTF)
 
 # obtain mukey map from SoilWeb Web Coverage Service (800m resolution SSURGO derived)
 res <- mukey.wcs(aoi = list(aoi = c(-114.16, 47.65,-114.08, 47.68), crs = 'EPSG:4326'))
-#> Loading required namespace: sf
 ```
 
-``` r
+    ## Loading required namespace: sf
 
+``` r
 # request input data from SDA
 varnames <- c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")
 resprop <- get_SDA_property(property = varnames,
                             method = "Dominant Component (numeric)",
                             mukeys = unique(values(res$mukey)))
+```
+
+    ## single result set, returning a data.frame
 
+``` r
 # keep only those where we have a complete set of 4 parameters (sand, silt, clay, bulk density; model code #3)
 soildata <- resprop[complete.cases(resprop), c("mukey", varnames)]
 
 # run Rosetta on the mapunit-level aggregate data
 system.time(resrose <- run_rosetta(soildata[,varnames]))
-#>    user  system elapsed 
-#>    0.03    0.00    0.06
 ```
 
-``` r
+    ##    user  system elapsed 
+    ##   0.022   0.005   0.027
 
+``` r
 # transfer mukey to result
 resprop$mukey <- as.numeric(resprop$mukey)
 resrose$mukey <- as.numeric(soildata$mukey)
@@ -282,7 +347,7 @@ res2 <- catalyze(res)
 plot(res2, "log10_Ksat_mean")
 ```
 
-<img src="man/figures/README-unnamed-chunk-10-1.png" width="100%" />
+![](README_files/figure-gfm/unnamed-chunk-11-1.png)<!-- -->
 
 ### *SpatRaster* (terra) Input Example
 
@@ -306,17 +371,17 @@ res3 <- rast(list(
 
 # SpatRaster to data.frame interface (one call on all cells)
 system.time(test2 <- run_rosetta(res3))
-#>    user  system elapsed 
-#>    6.20    0.55   14.42
 ```
 
-``` r
+    ##    user  system elapsed 
+    ##   5.312   8.671  13.721
 
+``` r
 # make a plot of the predicted Ksat (identical to mukey-based results)
 plot(test2, "log10_Ksat_mean")
 ```
 
-<img src="man/figures/README-unnamed-chunk-11-1.png" width="100%" />
+![](README_files/figure-gfm/unnamed-chunk-12-1.png)<!-- -->
 
 You will notice the results for Ksat distribution are identical since
 the same input values were used, but the latter approach took longer to
@@ -340,126 +405,88 @@ my_rosetta <- Rosetta(rosetta_version = 3, model_code = 3)
 
 ``` r
 predict(my_rosetta, list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
-#> [[1]]
-#>            [,1]     [,2]      [,3]      [,4]      [,5]
-#> [1,] 0.11535773 0.417912 -2.067139 0.1120102 0.8325407
-#> [2,] 0.09130753 0.485032 -2.022388 0.1510716 1.9060148
-#> 
-#> [[2]]
-#>            [,1]        [,2]       [,3]       [,4]       [,5]
-#> [1,] 0.01335011 0.009377977 0.08251142 0.01323413 0.09245277
-#> [2,] 0.01277141 0.013062171 0.10020312 0.01763982 0.14163567
 ```
 
+    ## $mean
+    ##            [,1]     [,2]      [,3]      [,4]      [,5]
+    ## [1,] 0.11535773 0.417912 -2.067139 0.1120102 0.8325407
+    ## [2,] 0.09130753 0.485032 -2.022388 0.1510716 1.9060148
+    ## 
+    ## $stdev
+    ##            [,1]        [,2]       [,3]       [,4]       [,5]
+    ## [1,] 0.01335011 0.009377977 0.08251142 0.01323413 0.09245277
+    ## [2,] 0.01277141 0.013062171 0.10020312 0.01763982 0.14163567
+
 ### Extended *Rosetta* Predictions, Parameter Distributions and Summary Statistics after Zhang & Schaap (2017) with `ann_predict()`
 
+`ann_predict()` is deprecated and redirects to `predict()`, as the
+underlying bootstrap data is now returned by `predict()` and summarized
+by R.
+
 ``` r
 ann_predict(my_rosetta, list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
-#> ann_predict() is defined for objects with class Rosetta; see `Rosetta()` to create a new instance
-#> $var_names
-#> $var_names[[1]]
-#> b'theta_r'
-#> 
-#> $var_names[[2]]
-#> b'theta_s'
-#> 
-#> $var_names[[3]]
-#> b'alpha'
-#> 
-#> $var_names[[4]]
-#> b'npar'
-#> 
-#> $var_names[[5]]
-#> b'ks'
-#> 
-#> 
-#> $sum_res_mean
-#>            [,1]        [,2]
-#> [1,]  0.1153577  0.09130753
-#> [2,]  0.4179120  0.48503196
-#> [3,] -2.0671385 -2.02238809
-#> [4,]  0.1120102  0.15107161
-#> [5,]  0.8325407  1.90601478
-#> 
-#> $sum_res_std
-#>             [,1]       [,2]
-#> [1,] 0.013350113 0.01277141
-#> [2,] 0.009377977 0.01306217
-#> [3,] 0.082511421 0.10020312
-#> [4,] 0.013234131 0.01763982
-#> [5,] 0.092452769 0.14163567
-#> 
-#> $sum_res_cov
-#> , , 1
-#> 
-#>              [,1]         [,2]          [,3]          [,4]         [,5]
-#> [1,] 1.782255e-04 3.719987e-05  5.300747e-05  2.278859e-05 5.965812e-05
-#> [2,] 3.719987e-05 8.794645e-05  2.570007e-04  3.632187e-06 2.662337e-04
-#> [3,] 5.300747e-05 2.570007e-04  6.808135e-03 -3.601399e-04 1.347519e-03
-#> [4,] 2.278859e-05 3.632187e-06 -3.601399e-04  1.751422e-04 8.973559e-05
-#> [5,] 5.965812e-05 2.662337e-04  1.347519e-03  8.973559e-05 8.547515e-03
-#> 
-#> , , 2
-#> 
-#>               [,1]          [,2]          [,3]          [,4]          [,5]
-#> [1,]  1.631088e-04  2.342660e-05  2.087559e-05 -9.477980e-06 -0.0003382175
-#> [2,]  2.342660e-05  1.706203e-04  3.156961e-04 -4.577936e-05  0.0005664825
-#> [3,]  2.087559e-05  3.156961e-04  1.004067e-02 -7.685234e-04  0.0012668368
-#> [4,] -9.477980e-06 -4.577936e-05 -7.685234e-04  3.111634e-04  0.0001828960
-#> [5,] -3.382175e-04  5.664825e-04  1.266837e-03  1.828960e-04  0.0200606627
-#> 
-#> 
-#> $sum_res_skew
-#>             [,1]        [,2]
-#> [1,] -4.52570431 -2.33577302
-#> [2,] -0.01729594 -0.22228088
-#> [3,] -0.18215435 -0.25961263
-#> [4,] -0.15345973  0.04145318
-#> [5,] -0.20386127 -0.35765020
-#> 
-#> $sum_res_kurt
-#>            [,1]        [,2]
-#> [1,] 36.6873026 16.85229297
-#> [2,]  0.5975976  0.02966211
-#> [3,]  0.3016512  0.11663125
-#> [4,]  0.1342581  0.26557316
-#> [5,]  0.3127817  0.54242091
-#> 
-#> $sum_res_bool
-#>      [,1] [,2]
-#> [1,] TRUE TRUE
-#> [2,] TRUE TRUE
-#> [3,] TRUE TRUE
-#> [4,] TRUE TRUE
-#> [5,] TRUE TRUE
-#> 
-#> $nsamp
-#> [1] 2
-#> 
-#> $nout
-#> [1] 5
-#> 
-#> $nin
-#> [1] 4
 ```
 
+    ## Warning in ann_predict.Rosetta(my_rosetta, list(c(30, 30, 40, 1.5), c(55, :
+    ## ann_predict() is deprecated in rosetta-soil >= 0.3.0. Use predict() instead.
+
+    ## $mean
+    ##            [,1]     [,2]      [,3]      [,4]      [,5]
+    ## [1,] 0.11535773 0.417912 -2.067139 0.1120102 0.8325407
+    ## [2,] 0.09130753 0.485032 -2.022388 0.1510716 1.9060148
+    ## 
+    ## $stdev
+    ##            [,1]        [,2]       [,3]       [,4]       [,5]
+    ## [1,] 0.01335011 0.009377977 0.08251142 0.01323413 0.09245277
+    ## [2,] 0.01277141 0.013062171 0.10020312 0.01763982 0.14163567
+
+## New Features in `rosetta-soil` 0.3.0
+
+### `rosesoil()`
+
+`rosesoil()` is a new R wrapper for the upstream `rosesoil()` function,
+which returns a structured result including all model metadata.
+
+``` r
+rosesoil(list(c(33, 33, 34, 1.5)))
+```
+
+    ##   sand silt clay rhob th33 th1500 version estimate_type code       thr
+    ## 1   33   33   34  1.5   NA     NA       3        linear    3 0.1076822
+    ##         ths      alpha     npar     ksat       k0       lpar    thr_std
+    ## 1 0.4060782 0.00817606 1.320054 7.141391 1.281258 -0.8606394 0.01221764
+    ##       ths_std   alpha_std   npar_std ksat_std    k0_std lpar_std
+    ## 1 0.008092908 0.001314515 0.03752358 1.415595 0.7639638 1.481137
+
+### `UnsaturatedK()`
+
+`UnsaturatedK` provides a way to predict unsaturated hydraulic
+conductivity parameters `K0` and `lpar` from retention parameters.
+
+``` r
+uk <- UnsaturatedK()
+predict(uk, list(c(0.12, 0.42, 0.008, 1.29)))
+```
+
+    ##   log10_K0_mean lpar_mean log10_K0_sd  lpar_sd
+    ## 1   -0.04057941 -1.033027   0.2309726 1.662214
+
 ## Selected References
 
 Three versions of the ROSETTA model are available, selected using
-`rosetta_version` argument.
+`rosetta_version` argument:
 
-- `rosetta_version` 1 - Schaap, M.G., F.J. Leij, and M.Th. van
+- `rosetta_version` 1: Schaap, M.G., F.J. Leij, and M.Th. van
   Genuchten. 2001. ROSETTA: a computer program for estimating soil
   hydraulic parameters with hierarchical pedotransfer functions. Journal
   of Hydrology 251(3-4): 163-176. doi: 10.1016/S0022-1694(01)00466-8.
-
-- `rosetta_version` 2 - Schaap, M.G., A. Nemes, and M.T. van
+- `rosetta_version` 2: Schaap, M.G., A. Nemes, and M.T. van
   Genuchten. 2004. Comparison of Models for Indirect Estimation of Water
   Retention and Available Water in Surface Soils. Vadose Zone Journal
   3(4): 1455-1463. doi: 10.2136/vzj2004.1455.
-
-- `rosetta_version` 3 - Zhang, Y., and M.G. Schaap. 2017. Weighted
+- `rosetta_version` 3: Zhang, Y., and M.G. Schaap. 2017. Weighted
   recalibration of the Rosetta pedotransfer model with improved
   estimates of hydraulic parameter distributions and summary statistics
   (Rosetta3). Journal of Hydrology 547: 39-53. doi:
-  10.1016/j.jhydrol.2017.01.004.
+  10.1016/j.jhydrol.2017.01.004. Version 3 includes predictions for
+  unsaturated conductivity parameters `K0` and `lpar`.
diff --git a/man/SoilDataFromArray.Rd b/man/SoilDataFromArray.Rd
index 9ebc52e..b974647 100644
--- a/man/SoilDataFromArray.Rd
+++ b/man/SoilDataFromArray.Rd
@@ -16,7 +16,7 @@ SoilDataFromArray(x)
 an object reference to a Rosetta \emph{SoilData} Python object constructed from \code{x}
 }
 \description{
-\code{SoilDataFromArray}: convert a list of numeric vectors containing soil properties to a \code{rosetta.rosetta.SoilData} class
+\code{SoilDataFromArray}: convert a list of numeric vectors containing soil properties to a \code{rosetta.rosetta.SoilData} class. In \code{rosetta-soil} >= 0.3, direct list input is preferred.
 
 \verb{py_to_r(<rosetta.rosetta.SoilData>)}: Wrapper S3 method for SoilData objects to prevent automatic conversion of SoilData (subclass of \code{"python.builtin.list"}) to an R \code{"list"}
 }
diff --git a/man/UnsaturatedK.Rd b/man/UnsaturatedK.Rd
new file mode 100644
index 0000000..619e720
--- /dev/null
+++ b/man/UnsaturatedK.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Class-Rosetta.R
+\name{UnsaturatedK}
+\alias{UnsaturatedK}
+\title{Make an UnsaturatedK object instance}
+\usage{
+UnsaturatedK()
+}
+\value{
+an instance of the \code{UnsaturatedK} class.
+}
+\description{
+\code{UnsaturatedK}: Create an instance of the \code{UnsaturatedK} class from \code{rosetta-soil} >= 0.3. This class is used to predict \code{K0} and \code{L} from retention parameters.
+}
diff --git a/man/figures/README-chunk-10-1.png b/man/figures/README-chunk-10-1.png
new file mode 100644
index 0000000..e6e8f43
Binary files /dev/null and b/man/figures/README-chunk-10-1.png differ
diff --git a/man/figures/README-chunk-11-1.png b/man/figures/README-chunk-11-1.png
new file mode 100644
index 0000000..42086f4
Binary files /dev/null and b/man/figures/README-chunk-11-1.png differ
diff --git a/man/figures/README-unnamed-chunk-10-1.png b/man/figures/README-unnamed-chunk-10-1.png
index e19f4bb..251474c 100644
Binary files a/man/figures/README-unnamed-chunk-10-1.png and b/man/figures/README-unnamed-chunk-10-1.png differ
diff --git a/man/figures/README-unnamed-chunk-11-1.png b/man/figures/README-unnamed-chunk-11-1.png
index 4e1cc4e..251474c 100644
Binary files a/man/figures/README-unnamed-chunk-11-1.png and b/man/figures/README-unnamed-chunk-11-1.png differ
diff --git a/man/figures/README-unnamed-chunk-12-1.png b/man/figures/README-unnamed-chunk-12-1.png
new file mode 100644
index 0000000..1af6490
Binary files /dev/null and b/man/figures/README-unnamed-chunk-12-1.png differ
diff --git a/man/predict.Rosetta.Rd b/man/predict.Rosetta.Rd
index f308043..168ef81 100644
--- a/man/predict.Rosetta.Rd
+++ b/man/predict.Rosetta.Rd
@@ -13,6 +13,12 @@
 
 \item{...}{not used}
 }
+\value{
+A list containing \code{mean} and \code{stdev} matrices (one row per sample).
+
+For \code{rosetta-soil} >= 0.3, the columns are: \code{theta_r}, \code{theta_s}, \code{alpha}, \code{npar}, \code{ksat}.
+Note that these parameters are in the scale produced by the underlying model (often log10 for alpha, npar, and ksat).
+}
 \description{
 Predict Rosetta Parameter Values and Standard Deviations from a \emph{Rosetta} instance
 }
diff --git a/man/predict.UnsaturatedK.Rd b/man/predict.UnsaturatedK.Rd
new file mode 100644
index 0000000..e2d5b10
--- /dev/null
+++ b/man/predict.UnsaturatedK.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Class-Rosetta.R
+\name{predict.UnsaturatedK}
+\alias{predict.UnsaturatedK}
+\title{Predict K0 and L from retention parameters}
+\usage{
+\method{predict}{UnsaturatedK}(object, retc_params, ...)
+}
+\arguments{
+\item{object}{\emph{UnsaturatedK} object}
+
+\item{retc_params}{A list or matrix of retention parameters (theta_r, theta_s, alpha, npar)}
+
+\item{...}{not used}
+}
+\value{
+a \code{data.frame} with \code{log10_K0_mean}, \code{lpar_mean}, \code{log10_K0_sd}, \code{lpar_sd}
+}
+\description{
+Predict K0 and L from retention parameters
+}
diff --git a/man/rosesoil.Rd b/man/rosesoil.Rd
new file mode 100644
index 0000000..a8e135f
--- /dev/null
+++ b/man/rosesoil.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rosesoil.R
+\name{rosesoil}
+\alias{rosesoil}
+\title{Run rosesoil() from rosetta-soil >= 0.3.0}
+\usage{
+rosesoil(soildata, rosetta_version = 3, estimate_type = "arith", vars = NULL)
+}
+\arguments{
+\item{soildata}{A list of numeric vectors or a data.frame (3-6 columns: sand, silt, clay, optionally bulk density, th33, and th1500)}
+
+\item{rosetta_version}{integer, 1-3. Default: 3}
+
+\item{estimate_type}{\emph{character}. One of \code{"arith"} (default), \code{"log"}, or \code{"geo"}. Only used if \code{rosetta-soil} >= 0.3.1. \code{"log"} returns parameters on a logarithmic (log10) scale for \code{alpha}, \code{npar}, \code{ksat}, and \code{k0}. \code{"geo"} returns the geometric mean of bootstrap estimates (exponent of the mean of log-transformed values). This is often preferred for parameters that vary by orders of magnitude, such as \code{alpha} and \code{ksat}.}
+
+\item{vars}{optional column name mapping (same as run_rosetta)}
+}
+\value{
+a data.frame with all RosettaResult fields
+}
+\description{
+Run rosesoil() from rosetta-soil >= 0.3.0
+}
diff --git a/man/rosettaPTF-package.Rd b/man/rosettaPTF-package.Rd
index acad19a..dce9b90 100644
--- a/man/rosettaPTF-package.Rd
+++ b/man/rosettaPTF-package.Rd
@@ -6,7 +6,7 @@
 \alias{rosettaPTF-package}
 \title{rosettaPTF: R Frontend for Rosetta Pedotransfer Functions}
 \description{
-Access Python rosetta-soil pedotransfer functions in an R environment. Rosetta is a neural network-based model for predicting unsaturated soil hydraulic parameters from basic soil characterization data. The model predicts parameters for the van Genuchten unsaturated soil hydraulic properties model, using sand, silt, and clay, bulk density and water content. The codebase is now maintained by Dr. Todd Skaggs and other U.S. Department of Agriculture employees. This R package is intended to provide for use cases that involve many thousands of calls to the pedotransfer function. Less demanding use cases are encouraged to use the web interface or API endpoint. There are additional wrappers of the API endpoints provided by the soilDB R package `ROSETTA()` method.
+Access the rosetta-soil Python pedotransfer functions from R. Rosetta is a neural network-based model for predicting unsaturated soil hydraulic parameters from basic soil characterization data (sand, silt, clay, bulk density, and water content). Predictions are made for the van Genuchten unsaturated hydraulic properties model, with uncertainty quantification via bootstrap ensemble. Designed for efficient batch processing of large datasets through vectorized computation and optional parallel processing.
 }
 \seealso{
 Useful links:
diff --git a/man/rosetta_pkg_version.Rd b/man/rosetta_pkg_version.Rd
new file mode 100644
index 0000000..51d1aa9
--- /dev/null
+++ b/man/rosetta_pkg_version.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rosetta_utils.R
+\name{rosetta_pkg_version}
+\alias{rosetta_pkg_version}
+\title{Get rosetta-soil Python package version}
+\usage{
+rosetta_pkg_version()
+}
+\value{
+\code{package_version} object
+}
+\description{
+Get rosetta-soil Python package version
+}
+\keyword{internal}
diff --git a/man/run_rosetta.Rd b/man/run_rosetta.Rd
index 4feba8b..0088aa2 100644
--- a/man/run_rosetta.Rd
+++ b/man/run_rosetta.Rd
@@ -10,16 +10,35 @@
 \alias{run_rosetta.SpatRaster}
 \title{Run \code{rosetta()} method from Python module}
 \usage{
-\method{run_rosetta}{default}(soildata, vars = NULL, rosetta_version = 3, ...)
+\method{run_rosetta}{default}(
+  soildata,
+  vars = NULL,
+  rosetta_version = 3,
+  estimate_type = "log",
+  ...
+)
 
-\method{run_rosetta}{data.frame}(soildata, vars = NULL, rosetta_version = 3, ...)
+\method{run_rosetta}{data.frame}(
+  soildata,
+  vars = NULL,
+  rosetta_version = 3,
+  estimate_type = "log",
+  ...
+)
 
-\method{run_rosetta}{matrix}(soildata, vars = NULL, rosetta_version = 3, ...)
+\method{run_rosetta}{matrix}(
+  soildata,
+  vars = NULL,
+  rosetta_version = 3,
+  estimate_type = "log",
+  ...
+)
 
 \method{run_rosetta}{RasterStack}(
   soildata,
   vars = NULL,
   rosetta_version = 3,
+  estimate_type = "log",
   cores = 1,
   core_thresh = 20000L,
   file = paste0(tempfile(), ".tif"),
@@ -31,6 +50,7 @@
   soildata,
   vars = NULL,
   rosetta_version = 3,
+  estimate_type = "log",
   cores = 1,
   core_thresh = 20000L,
   file = paste0(tempfile(), ".tif"),
@@ -42,6 +62,7 @@
   soildata,
   vars = NULL,
   rosetta_version = 3,
+  estimate_type = "log",
   cores = 1,
   core_thresh = 20000L,
   file = paste0(tempfile(), ".tif"),
@@ -56,6 +77,8 @@
 
 \item{rosetta_version}{Default: 3}
 
+\item{estimate_type}{\emph{character}. One of \code{"log"} (default), \code{"arith"}, or \code{"geo"}. Only used if \code{rosetta-soil} >= 0.3.1. Default \code{"log"} preserves logarithmic (log10) scale for \code{alpha}, \code{npar}, and \code{Ksat}. \code{"geo"} returns the geometric mean of bootstrap estimates (exponent of the mean of log-transformed values). This is often preferred for parameters that vary by orders of magnitude, such as \code{alpha} and \code{Ksat}.}
+
 \item{...}{additional arguments not used}
 
 \item{cores}{number of cores; used only for processing \emph{SpatRaster} or \emph{Raster*} input}
@@ -69,13 +92,15 @@
 \item{overwrite}{logical; overwrite \code{file}? passed to \code{terra::writeStart()}; defaults to \code{TRUE} if needed}
 }
 \value{
-A \emph{data.frame} containing \code{mean} and \code{stdev} for following five columns (parameters for van Genuchten-Mualem equation)
+A \emph{data.frame} containing \code{mean} and \code{stdev} for the following columns (parameters for van Genuchten-Mualem equation)
 \itemize{
 \item \code{"theta_r"}, residual water content
 \item \code{"theta_s"}, saturated water content
-\item \code{"log10(alpha)"}, 'alpha' shape parameter, log10(1/cm)
-\item \code{"log10(npar)"}, 'n' shape parameter
-\item \code{"log10(Ksat)"}, saturated hydraulic conductivity, log10(cm/day)
+\item \code{"alpha"}, 'alpha' shape parameter (1/cm). Logarithmic (log10) scale if \code{estimate_type="log"} (default); Geometric mean if \code{estimate_type="geo"}.
+\item \code{"npar"}, 'n' shape parameter. Logarithmic (log10) scale if \code{estimate_type="log"} (default); Geometric mean if \code{estimate_type="geo"}.
+\item \code{"Ksat"}, saturated hydraulic conductivity (cm/day). Logarithmic (log10) scale if \code{estimate_type="log"} (default); Geometric mean if \code{estimate_type="geo"}.
+\item \code{"K0"}, unsaturated hydraulic conductivity (cm/day). Only if \code{rosetta-soil} >= 0.3.1. Logarithmic (log10) scale if \code{estimate_type="log"} (default); Geometric mean if \code{estimate_type="geo"}.
+\item \code{"lpar"}, unsaturated hydraulic conductivity exponent. Only if \code{rosetta-soil} >= 0.3.1.
 }
 
 If the sum of sand, silt, and clay is not 100\%, the parameter value estimates will be \code{NaN}.
@@ -83,3 +108,9 @@ If the sum of sand, silt, and clay is not 100\%, the parameter value estimates w
 \description{
 Run \code{rosetta()} method from Python module
 }
+\details{
+\subsection{Performance Note}{
+
+Use \code{cores > 1} with \code{SpatRaster} or \verb{Raster*} inputs to parallelize processing of cells across multiple cores.
+}
+}
diff --git a/tests/testthat/test-predict-Rosetta.R b/tests/testthat/test-predict-Rosetta.R
index cdb2835..302a741 100644
--- a/tests/testthat/test-predict-Rosetta.R
+++ b/tests/testthat/test-predict-Rosetta.R
@@ -4,7 +4,20 @@ test_that("prediction with Rosetta class works", {
   skip_if_not(py_module_available("rosetta"))
 
   one <- predict(Rosetta(), list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
-  two <- ann_predict(Rosetta(), list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1)))
+  expect_warning({ two <- ann_predict(Rosetta(), list(c(30, 30, 40, 1.5), c(55, 25, 20, 1.1))) })
 
   expect_true(inherits(one, 'list') && inherits(two, 'list'))
+  expect_true("mean" %in% names(one))
+  expect_true("stdev" %in% names(one))
+})
+
+test_that("UnsaturatedK works", {
+  skip_if_not(py_module_available("numpy"))
+  skip_if_not(py_module_available("rosetta"))
+  skip_if(rosetta_pkg_version() < package_version("0.3.0"))
+
+  uk <- UnsaturatedK()
+  res <- predict(uk, list(c(0.12, 0.42, 0.008, 1.29)))
+  expect_true(inherits(res, "data.frame"))
+  expect_true("log10_K0_mean" %in% colnames(res))
 })
diff --git a/tests/testthat/test-rosesoil.R b/tests/testthat/test-rosesoil.R
new file mode 100644
index 0000000..bbeac78
--- /dev/null
+++ b/tests/testthat/test-rosesoil.R
@@ -0,0 +1,21 @@
+test_that("rosesoil() works", {
+  skip_if_not(py_module_available("rosetta"))
+  skip_if(rosetta_pkg_version() < package_version("0.3.0"))
+  
+  res <- rosesoil(list(c(30, 30, 40, 1.5)))
+  expect_true(inherits(res, "data.frame"))
+  expect_true("thr" %in% colnames(res))
+  expect_true("ths" %in% colnames(res))
+  expect_true("k0" %in% colnames(res))
+  expect_true("lpar" %in% colnames(res))
+})
+
+test_that("rosesoil() with data.frame and vars works", {
+  skip_if_not(py_module_available("rosetta"))
+  skip_if(rosetta_pkg_version() < package_version("0.3.0"))
+  
+  df <- data.frame(S = 30, Si = 30, C = 40, BD = 1.5)
+  res <- rosesoil(df, vars = c("S", "Si", "C", "BD"))
+  expect_true(inherits(res, "data.frame"))
+  expect_true("thr" %in% colnames(res))
+})
diff --git a/tests/testthat/test-rosetta.R b/tests/testthat/test-rosetta.R
index 0cdcd3b..3698541 100644
--- a/tests/testthat/test-rosetta.R
+++ b/tests/testthat/test-rosetta.R
@@ -1,32 +1,56 @@
-test_that("run_rosetta() works", {
+test_that("run_rosetta() with sample data", {
   skip_if_not(py_module_available("numpy"))
   skip_if_not(py_module_available("rosetta"))
 
-  res <- run_rosetta(list(c(30, 30, 40, 1.5), c(55, 25, 20),  c(55, 25, 20, 1.1)),
-                     rosetta_version = 3)
+  data("MUKEY_PROP")
+  varnames <- c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")
+  res <- run_rosetta(MUKEY_PROP[1:10, varnames], rosetta_version = 3)
+
   expect_true(inherits(res, 'data.frame'))
+  if (rosetta_pkg_version() >= package_version("0.3.0")) {
+    expect_true(ncol(res) >= 15)
+    expect_true("log10_K0_mean" %in% colnames(res))
+  } else {
+    expect_true(ncol(res) == 12)
+  }
 })
 
-test_that("data.frame interface", {
+test_that("estimate_type argument with sample data", {
+  skip_if_not(py_module_available("rosetta"))
+  skip_if(rosetta_pkg_version() < package_version("0.3.0"))
+
+  data("MUKEY_PROP")
+  varnames <- c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")
+
+  res_log <- run_rosetta(MUKEY_PROP[1:5, varnames], estimate_type = "log")
+  expect_true("log10_Ksat_mean" %in% colnames(res_log))
 
+  res_lin <- run_rosetta(MUKEY_PROP[1:5, varnames], estimate_type = "arith")
+  expect_true("ksat_mean" %in% colnames(res_lin))
+  expect_false("log10_Ksat_mean" %in% colnames(res_lin))
+
+  res_geo <- run_rosetta(MUKEY_PROP[1:5, varnames], estimate_type = "geo")
+  expect_true("ksat_mean" %in% colnames(res_geo))
+})
+
+test_that("data.frame interface", {
   skip_if_not(py_module_available("numpy"))
   skip_if_not(py_module_available("rosetta"))
 
-  # data.frame interface: using default column order
-  expect_true(inherits(run_rosetta(data.frame(
-    a = 20,
-    b = 60,
-    c = 20,
-    d = c(NA, 1.5)
-  )), 'data.frame'))
+  # Default column order with sample data
+  data("MUKEY_PROP")
+  varnames <- c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")
+  res1 <- run_rosetta(MUKEY_PROP[1:5, varnames])
+  expect_true(inherits(res1, 'data.frame'))
 
-  # data.frame interface: using custom column names/order
-  expect_true(inherits(run_rosetta(data.frame(
+  # Custom column order with synthetic data
+  res2 <- run_rosetta(data.frame(
     d = c(NA, 1.5),
     b = 60,
     a = 20,
     c = 20
-  ), vars = letters[1:4]), 'data.frame'))
+  ), vars = letters[1:4])
+  expect_true(inherits(res2, 'data.frame'))
 })
 
 test_that("run on SSURGO data", {
@@ -34,8 +58,8 @@ test_that("run on SSURGO data", {
   skip_if_not(py_module_available("rosetta"))
 
   data("MUKEY_WCS", package = "rosettaPTF")
-  res <- terra::rast(MUKEY_WCS, crs = "EPSG:6350")
-  terra::ext(res) <- c(-114.16, 47.65, -114.08, 47.68)
+  res <- terra::rast(MUKEY_WCS, crs = "EPSG:5070")
+  terra::ext(res) <- c(-1365495, -1358925, 2869245, 2873655)
   names(res) <- "mukey"
 
   mukeys <- as.numeric(terra::values(res[[1]]))
diff --git a/vignettes/performance-raster.Rmd b/vignettes/performance-raster.Rmd
new file mode 100644
index 0000000..534483f
--- /dev/null
+++ b/vignettes/performance-raster.Rmd
@@ -0,0 +1,155 @@
+---
+title: "High-Throughput Raster Processing"
+knit: litedown:::knit
+vignette: >
+  %\VignetteIndexEntry{High-Throughput Raster Processing}
+  %\VignetteEngine{litedown::vignette}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+library(rosettaPTF)
+library(terra)
+
+EVAL <- rosettaPTF::rosetta_module_available()
+
+litedown::reactor(
+  eval = EVAL,
+  collapse = TRUE,
+  fig.width = 8,
+  fig.align = 'center'
+)
+```
+
+## Introduction
+
+`{rosettaPTF}` is designed for efficient batch processing of soil hydraulic parameters. The package supports multiple input formats and offers options for scaling analysis to large datasets.
+
+The core implementation uses the `rosetta-soil` Python module, which provides vectorized computation over multiple samples. This vignette demonstrates key parameters for controlling performance and handling large-scale analyses.
+
+## Processing Point Data
+
+The `run_rosetta()` function accepts point data as a `data.frame`, where each row represents a single observation. Let's use the sample soil property dataset included with the package.
+
+```{r point_data}
+data("MUKEY_PROP")
+
+# View the structure of the sample data
+str(MUKEY_PROP)
+
+# Run rosetta on the sample property data
+system.time({
+  res_points <- run_rosetta(MUKEY_PROP[, c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")])
+})
+
+head(res_points)
+```
+
+## Processing Continuously Varying Raster Data
+
+Spatial soil property predictions are commonly stored as raster grids. The `run_rosetta()` function can process `SpatRaster` objects directly, computing predictions for every cell in the grid.
+
+### Creating a Raster Stack from Sample Data
+
+We'll use the sample spatial dataset (MUKEY_WCS) to create a continuous raster surface by interpolating soil properties.
+
+```{r raster_setup}
+# Convert MUKEY_WCS matrix to SpatRaster
+data("MUKEY_WCS", package = "rosettaPTF")
+data("MUKEY_PROP", package = "rosettaPTF")
+
+r_template <- terra::rast(MUKEY_WCS, crs = "EPSG:5070")
+terra::ext(r_template) <- c(-1365495, -1358925, 2869245, 2873655)
+names(r_template) <- "mukey"
+
+levels(r_template) <- MUKEY_PROP[, c("mukey",
+  "sandtotal_r", "silttotal_r", "claytotal_r",
+  "dbthirdbar_r")]
+
+r_input <- terra::catalyze(r_template)
+plot(r_input)
+```
+
+### Running Rosetta on Raster Data
+
+Pass the `SpatRaster` object to `run_rosetta()`. The output is a multi-layer raster containing mean and standard deviation for each predicted parameter.
+
+```{r raster_processing}
+# Process the raster stack
+system.time({
+  r_output <- run_rosetta(r_input)
+})
+
+# Inspect the layers
+names(r_output)
+
+# Plot predicted Ksat (log10 cm/day)
+plot(r_output[["log10_Ksat_mean"]], main = "Predicted Ksat")
+```
+
+## Scaling to Large Datasets
+
+### Parallel Processing with Multiple Cores
+
+For large rasters, the `cores` argument enables block-wise processing across multiple CPU cores. This parameter controls how the raster is divided and processed in parallel.
+
+```{r parallel_raster, eval=FALSE}
+# Divide the raster into blocks and process each block on separate cores
+r_output_parallel <- run_rosetta(r_input, cores = 2)
+```
+
+Note that with small rasters, as in this example, parallel processing may be significantly slower than sequential processing.
+
+### Key Parameters for Scaling
+
+- **`cores`**: Number of CPU cores to use for parallel processing. Set to 1 (default) for sequential processing, or 2+ for parallel block-wise processing. Useful for rasters that are memory-intensive or computationally demanding.
+
+- **Input format**: `SpatRaster` objects are preferred for spatial workflows. The function handles memory-efficient extraction and output reconstruction automatically.
+
+### Batch Processing Workflows
+
+For extremely large regions or high-resolution grids, consider processing tiles or regions sequentially:
+
+```{r batch_example, eval=FALSE}
+# Example: process raster in regional tiles
+tiles <- terra::getTileExtents(r_input, 125)
+
+results <- terra::merge(terra::sprc(apply(tiles, 1, function(x) {
+  terra::window(r_input) <- x
+  run_rosetta(r_input)
+})))
+```
+
+## Advanced Options
+
+### Controlling Parameter Estimation Scale
+
+By default, parameters like `alpha`, `npar`, and `Ksat` are returned on a logarithmic (log10) scale. Alternative scales are available via the `estimate_type` argument:
+
+```{r estimate_scales}
+# Linear scale estimates
+res_linear <- run_rosetta(MUKEY_PROP[, c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")],
+                          estimate_type = "arith")
+
+# Geometric mean (recommended for log-transformed parameters)
+res_geo <- run_rosetta(MUKEY_PROP[, c("sandtotal_r", "silttotal_r", "claytotal_r", "dbthirdbar_r")],
+                       estimate_type = "geo")
+```
+
+### Bootstrap Ensemble
+
+By default, predictions include the full 1,000-member bootstrap ensemble for uncertainty quantification. The output includes both mean and standard deviation for each parameter.
+
+```{r uncertainty}
+# The output includes uncertainty estimates
+head(res_points[, c("log10_Ksat_mean", "log10_Ksat_sd")])
+```
+
+## Summary
+
+Key considerations for high-throughput analysis:
+
+1. Use `SpatRaster` objects as input for spatial workflows; the function handles data conversion automatically.
+2. Set `cores > 1` to enable parallel block-wise processing for large rasters.
+3. Choose `estimate_type` based on your application requirements.
+4. For extremely large regions, consider tiling or regional batch processing workflows.