R/watbal_basin_of.R

Defines functions watbal_basin_of

Documented in watbal_basin_of

#' watbal_basin_of.R
#'
#' Water balance for basin daily RHESSys output, as generated by outputfilters. To run, if your basin output is called "bd" type "bd=watbal(bd)".
#' This will add a number of fields to your basin output file, the last of which will be called "watbal".
#' This is your water balance error. You should see a minor numerical error here even if your water balances
#' (on the order of 10^-6).  If your watbal values are negative then water inputs are less than water outputs, and vice versa.
#' @param bd The basin daily outputs from rhessys, most easily retrieved via `readin_rhessys_output()`
#' @export

watbal_basin_of = function(bd) {

  # some error checks here
  req_cols = c("total_water_in", "evaporation", "evaporation_surf", "exfiltration_sat_zone",
  "exfiltration_unsat_zone", "transpiration_sat_zone", "transpiration_unsat_zone", "streamflow", "gw.Qout", "sat_deficit", "rz_storage", "unsat_storage",
	 "detention_store", "litter.rain_stored", "gw.storage","canopy_rain_stored", "canopy_snow_stored","snowpack.water_equivalent_depth")

  bd$evap = with(bd, evaporation+evaporation_surf+exfiltration_sat_zone+exfiltration_unsat_zone)
  bd$trans = with(bd, transpiration_unsat_zone+transpiration_sat_zone )
  bd$streamflow = with(bd, streamflow+gw.Qout)
  bd$canopy_store = with(bd, canopy_snow_stored+canopy_rain_stored )
  if (!is.data.frame(bd) || any(!req_cols %in% colnames(bd))) {
    cat("Input is either not a data frame or is missing the correct columns")
    return(NA)
  }

  # main fluxes
  bd$watbal_flux = with(bd, total_water_in - streamflow - trans - evap)

  # changes in stores
  bd$sd = with(bd, sat_deficit - rz_storage - unsat_storage)
  bd$sddiff = c(0, diff(bd$sd))
  bd$snodiff = c(0, diff(bd$snowpack.water_equivalent_depth))
  bd$detdiff = c(0, diff(bd$detention_store))
  bd$litdiff = c(0, diff(bd$litter.rain_stored))
  bd$candiff = c(0, diff(bd$canopy_store))
  bd$gwdiff = c(0, diff(bd$gw.storage))

  # fluxes minus stores
  bd$watbal = with(bd, watbal_flux + sddiff - snodiff - detdiff - litdiff - candiff - gwdiff)
  bd$watbal[1] = 0.0
  # max(bd$watbal)
  # summary(bd$watbal)
  # hist(bd$watbal)

  return(bd)

}
ryanrbart/rhessysR documentation built on March 29, 2024, 4:30 p.m.