BASICSimulationModel: proposed changes to cryic/frigid logic

[brownag]

• Newhall model gains support for additional spatially explicit inputs:

  • O horizon (presence/absence)
  • Saturation (presence/absence)
  • soilAirOffset (degrees C or degrees F)

• Cryic Soil Temperature Regime criteria now properly incorporate O horizon and saturation into thresholds

---

DRAFT PR for discussion and testing of changes to cryic/frigid determination in BASICSimulationModel

• it appears one of the cryic criteria was misinterpreted at some point, will try to translate it literally and see if results change
• as background (more to come on this) climate simulations run in the Sierra Nevada anecdotally skip from mesic directly to cryic. Need to evaluate which of the criteria are triggering this.
• May fuel some future proposals for clarifying how to apply cryic criteria 1a and 1b on p.31 in KST (12th edition)

[brownag]

library(jNSMR)
#> jNSMR (0.3.0) -- R interface to the classic Newhall Simulation Model
#> Added JAR file (newhall-1.6.5.jar) to Java class path.
library(soilDB)
library(terra)
library(rgeedim)

x <- vect(fetchSDA_spatial(c("CA067", "CA628", "CA731", "CA729", "CA732"),
                      by.col = "areasymbol",
                      geom.src = "sapolygon"))
#> Using 1 chunks...
#> Chunk #1 completed (n = 5; 4.2 secs)
#> Done in 5 secs; mean/chunk: 4.2 secs; mean/symbol: 1 secs.

# get PRISM climate and SSURGO-derived available water storage data
predictors <- c(newhall_prism_subset(x), newhall_issr800_subset(x))

# take focal window to smooth/remove NoData in AWC
predictors$awc <- focal(predictors$awc, 5, fun="mean", na.rm=TRUE)

# use Lang et al. (2023) canopy height <https://www.nature.com/articles/s41559-023-02206-6>
gd_initialize()

gd_image_from_id('users/nlang/ETH_GlobalCanopyHeight_2020_10m_v1') |>
  gd_download(region = as.polygons(x, ext = TRUE),
              bands = list("b1"),
              filename="canopy_height.tif",
              scale = 100)
#> [1] "canopy_height.tif"
h <- classify(rast("canopy_height.tif"), cbind(NA, 0))

# inspect
plot(h)

# assume if canopy is 2x minimum for tree height (4m) there is an O horizon
predictors$hasOHorizon <- project(h, predictors) > 8

plot(predictors$hasOHorizon)

# add elevation
gd_image_from_id('USGS/NED') |>
  gd_download(region = as.polygons(x, ext = TRUE),
              bands = list("elevation"),
              filename="elev.tif",
              scale = 100)
#> [1] "elev.tif"
e <- rast("elev.tif")
predictors$elev <- project(e, predictors)

# run batches of models on gridded data
res <- newhall_batch(crop(predictors, ext(predictors)), cores=8, core_thresh = 1e5)
#> Batch 1 of 2 (n=99660 on 8 cores) done in 1.3 mins
#> Batch 2 of 2 (n=61152 on 8 cores) done in 40 secs

y <- fetchSDA_spatial(SDA_spatialQuery(predictors, "areasymbol")$areasymbol,
                      by.col = "areasymbol",
                      geom.src = "sapolygon")
#> Using 6 chunks...
#> Chunk #1 completed (n = 10; 5.3 secs)
#> Chunk #2 completed (n = 10; 9.4 secs)
#> Chunk #3 completed (n = 10; 3 secs)
#> Chunk #4 completed (n = 10; 5.1 secs)
#> Chunk #5 completed (n = 10; 8.1 secs)
#> Chunk #6 completed (n = 4; 0.3 secs)
#> Done in 32.2 secs; mean/chunk: 5.2 secs; mean/symbol: 0.6 secs.
# inspect one modeled soil taxonomic property
plot(
  res$dryDaysAfterSummerSolstice,
  main = "dryDaysAfterSummerSolstice",
  col = grDevices::hcl.colors(50, "cividis")
)
plot(as.lines(vect(y)), add = TRUE, col = "darkorange")

# apply xeric threshold (dry >1.5 months after summer solstice)
plot(
  res$dryDaysAfterSummerSolstice > 45,
  main = "dryDaysAfterSummerSolstice",
  col = grDevices::hcl.colors(2, "cividis")
)
plot(as.lines(vect(y)), add = TRUE, col = "darkorange")

plot(res$temperatureRegime,
     col = grDevices::hcl.colors(7, "cividis"))
plot(as.lines(vect(y)), add = TRUE, col = "darkorange")