Water Balance Functions - Rwema25/AE-project GitHub Wiki

Overview of the Soil Water Balance Function Chain and Workflow

The soil water balance calculation is structured around a core function, calc_daily_watbal, which integrates several computation steps implemented through subfunctions (AWCPTF, soilcap_calc, calc_pet, eabyep_subfun, calc_watbal_series, and validate_watbal ) and internal calculations. The workflow starts with soil and weather input parameters and ends with a detailed daily water balance output.

The process begins with reading and processing both soil and weather data to prepare a clean and integrated dataset for the full model function described below. Soil horizon data are sourced from two datasets (ISDA and SoilGrids2.0), while weather data are obtained from the POWER and CHIRPS datasets. The following table summarizes the key variables used from each source. R code used

Soil Horizon Data		Weather Data
ISDA	SoilGrids2.0	POWER	CHIRPS
- sand.tot.psa: Sand content - silt.tot.psa: Silt content - clay.tot.psa: Clay content - oc: Soil organic carbon - db.od: Bulk density - ecec.f: Cation exchange capacity - ph.h2o: Soil pH in H2O	- sand: Sand content - silt: Silt content - clay: Clay content - soc: Soil organic carbon - bdod: Bulk density - cec: Cation exchange capacity - phh2o: Soil pH in H2O	- srad: Solar radiation - tmin: Minimum temperature - tmax: Maximum temperature - tmean: Mean temperature	- Rainfall

The following section outlines the subfunctions developed to support the main run_full_water_balance function, highlighting their specific roles in the workflow, detailing their input and output variables, and providing the corresponding R code used.

1. AWCPTF Function (Available Water Capacity Pedotransfer Function)

Description

Pedo-transfer function for available water capacity (AWC), based on Hodnett and Tomasella (2002). Calculates volumetric water content at various suctions (e.g., field capacity, wilting point), and from these values, derives available water capacity.

Input variable	Output variable	Data source
- SNDPPT - SLTPPT - CLYPPT - ORCDRC - BLD = 1400 - CEC - PHIHOX - h1 = -10 - h2 = -20 - h3 = -31.6 - pwp = -1585 - PTF.coef	- AWC at suction h1 - AWC at suction h2 - AWC at suction h3 (commonly used for field capacity) - Volumetric water content at wilting point (-1585 kPa), - Saturated volumetric water content (tetaS)	- SoilGrids2.0 - ISDA

Code

# Author: Peter Steward, Fabio Castro, Julian Rameriz
# Date: 2025-03-15

AWCPTF <- function(
    SNDPPT,
    SLTPPT,
    CLYPPT,
    ORCDRC,
    BLD = 1400,
    CEC,
    PHIHOX,
    h1 = -10,
    h2 = -20,
    h3 = -31.6,
    pwp = -1585,
    PTF.coef,
    fix.values = TRUE,
    print.coef = TRUE
){
  # If no coefficient table given, use defaults from Hodnett & Tomasella (2002).
  if (missing(PTF.coef)){
    PTF.coef <- data.frame(
      lnAlfa = c(-2.294,  0, -3.526,  0,    2.44,  0, -0.076, -11.331,  0.019,  0,     0,     0),
      lnN    = c(62.986,  0,  0,      -0.833,-0.529,0,  0,      0.593,   0,     0.007, -0.014, 0),
      tetaS  = c(81.799,  0,  0,       0.099, 0,    -31.42,0.018, 0.451,   0,     0,      0,     -5e-04),
      tetaR  = c(22.733, -0.164,0,     0,     0,     0,     0.235, -0.831,  0,     0.0018, 0,     0.0026)
    )
  }
  
  if (fix.values){
    sum.tex <- CLYPPT + SLTPPT + SNDPPT
    # Convert to percentages if needed
    CLYPPT <- CLYPPT / sum.tex * 100
    SLTPPT <- SLTPPT / sum.tex * 100
    SNDPPT <- SNDPPT / sum.tex * 100
    BLD[BLD < 100]  <- 100
    BLD[BLD > 2650] <- 2650
  }
  
  clm <- data.frame(
    SNDPPT,
    SLTPPT,
    CLYPPT,
    ORCDRC = ORCDRC / 10,
    BLD    = BLD * 0.001,
    CEC,
    PHIHOX,
    SLTPPT_sq = SLTPPT^2,
    CLYPPT_sq = CLYPPT^2,
    SNDxSLT   = SNDPPT * SLTPPT,
    SNDxCLY   = SNDPPT * CLYPPT
  )
  
  alfa <- apply(clm, 1, function(x) {
    exp((PTF.coef$lnAlfa[1] + sum(PTF.coef$lnAlfa[-1] * x)) / 100)
  })
  
  N <- apply(clm, 1, function(x) {
    exp((PTF.coef$lnN[1] + sum(PTF.coef$lnN[-1] * x)) / 100)
  })
  
  tetaS <- apply(clm, 1, function(x) {
    (PTF.coef$tetaS[1] + sum(PTF.coef$tetaS[-1] * x)) / 100
  })
  
  tetaR <- apply(clm, 1, function(x) {
    (PTF.coef$tetaR[1] + sum(PTF.coef$tetaR[-1] * x)) / 100
  })
  
  # Fix extremes
  tetaR[tetaR < 0]   <- 0
  tetaS[tetaS > 100] <- 100
  
  # Compute volumetric water content at the given suctions
  m       <- 1 - 1/N
  tetah1  <- tetaR + (tetaS - tetaR) / ((1 + (alfa * -1 * h1)^N))^m
  tetah2  <- tetaR + (tetaS - tetaR) / ((1 + (alfa * -1 * h2)^N))^m
  tetah3  <- tetaR + (tetaS - tetaR) / ((1 + (alfa * -1 * h3)^N))^m
  WWP     <- tetaR + (tetaS - tetaR) / ((1 + (alfa * -1 * pwp)^N))^m
  
  # If WWP is larger than one of the tetahs, clamp
  if (fix.values){
    sel <- which(WWP > tetah1 | WWP > tetah2 | WWP > tetah3)
    if (length(sel) > 0){
      WWP[sel] <- apply(
        data.frame(tetah1[sel], tetah2[sel], tetah3[sel]),
        1,
        min,
        na.rm=TRUE
      )
      warning(paste("Wilting point capacity is higher than h1/h2/h3 for",
                    length(sel), "points. Adjusted."))
    }
  }
  
  # Available Water Content for each suction
  AWCh1 <- tetah1 - WWP
  AWCh2 <- tetah2 - WWP
  AWCh3 <- tetah3 - WWP
  
  out <- data.frame(
    AWCh1 = signif(AWCh1, 3),
    AWCh2 = signif(AWCh2, 3),
    AWCh3 = signif(AWCh3, 3),
    WWP   = signif(WWP,   3),
    tetaS = signif(tetaS, 3)
  )
  
  if (print.coef){
    attr(out, "coef")      <- as.list(PTF.coef)
    attr(out, "PTF.names") <- list(
      variable = c("ai1","sand","silt","clay","oc","bd","cec","ph","silt^2",
                   "clay^2","sand*silt","sand*clay")
    )
  }
  
  return(out)
}

2. soilcap_calc Function

Description

Computes total water-holding capacity (mm) up to a given rooting depth, based on horizon-wise AWC values and lower-bound depths.

Input variable	Output variable	Data source
- x: A numeric vector of volumetric fractions (such as Available Water Content values for each soil horizon. e.g., AWCh3) - y: A numeric vector of horizon lower depth limits in centimeters - rdepth: Numeric value indicating the rooting depth in centimeters (default is 60 cm) - minval: Numeric minimum plausible rooting depth (in cm) - maxval: Numeric maximum plausible rooting depth (in cm)	- Total water holding capacity (mm) from surface to rdepth (soilcp)	- AWCPTF

Code

# Author: Peter Steward, Fabio Castro, Julian Rameriz
# Date: 2025-03-15

soilcap_calc <- function(x, y, rdepth=60, minval=45, maxval=100){
  if (length(x) != length(y)){
    stop("Length of x and y must be the same.")
  }
  # Force rdepth within [minval, maxval]
  rdepth <- max(c(rdepth, minval))
  rdepth <- min(c(rdepth, maxval))
  
  wc_df <- data.frame(depth=y, wc=x)
  
  # Interpolate if needed
  if (!rdepth %in% wc_df$depth){
    wc_df1 <- wc_df[wc_df$depth < rdepth, ]
    wc_df2 <- wc_df[wc_df$depth > rdepth, ]
    if (nrow(wc_df1) > 0 && nrow(wc_df2) > 0){
      y1 <- wc_df1$wc[nrow(wc_df1)]
      y2 <- wc_df2$wc[1]
      x1 <- wc_df1$depth[nrow(wc_df1)]
      x2 <- wc_df2$depth[1]
      ya <- (rdepth - x1)/(x2 - x1)*(y2 - y1) + y1
      wc_df <- rbind(
        wc_df1,
        data.frame(depth=rdepth, wc=ya),
        wc_df2
      )
    }
    # If the entire soil is shallower or deeper than rdepth,
    # the standard logic proceeds anyway.
  }
  
  wc_df <- wc_df[wc_df$depth <= rdepth, ]
  wc_df$soilthick <- wc_df$depth - c(0, wc_df$depth[-nrow(wc_df)])
  wc_df$soilcap   <- wc_df$soilthick * wc_df$wc
  
  # Convert from cm * fraction to mm
  soilcp <- sum(wc_df$soilcap) * 10
  return(soilcp)
}

3. calc_pet Function (Potential Evapotranspiration Function)

Description

Calculates potential evapotranspiration (PET, mm/day) using a Priestley-Taylor approach, which is a simplified model for estimating ET based on climatic variables.

Input variable	Output variable	Data source
- srad: Solar radiation (MJ/m2/day) - tmin: Minimum temperature (degree Celsius) - tmax: Maximum temperature (degree Celsius) - tmean: Mean temperature (degree Celsius)	- PET: Potential evapotranspiration (PET, mm/day)	- POWER

code

# Author: Peter Steward, Fabio Castro, Julian Rameriz
# Date: 2025-03-15

calc_pet <- function(
    srad,
    tmin,
    tmax,
    tmean=NULL,
    albedo=0.2,
    vpd_cte=0.7,
    psycho=62,
    rlat_ht=2.26e6,
    rho_w=997
){
  if (is.null(tmean)){
    tmean <- (tmin + tmax)/2
  }
  # Net radiation
  rn <- (1 - albedo)*srad
  
  # Slope of saturation vapor pressure curve
  a_eslope <- 611.2
  b_eslope <- 17.67
  c_eslope <- 243.5
  eslope <- a_eslope*b_eslope*c_eslope/(tmean + c_eslope)^2 * exp(b_eslope*tmean/(tmean + c_eslope))
  
  # Estimate VPD
  esat_min <- 0.61120*exp((17.67*tmin)/(tmin+243.5))
  esat_max <- 0.61120*exp((17.67*tmax)/(tmax+243.5))
  vpd      <- vpd_cte * (esat_max - esat_min)  # kPa
  
  # Priestley-Taylor factor
  pt_const <- 1.26
  vpd_ref  <- 1
  pt_fact  <- 1
  pt_coef  <- pt_fact * pt_const
  pt_coef2 <- 1 + (pt_coef - 1) * (vpd / vpd_ref)
  
  # Convert flux to mm water
  et_max <- (pt_coef2*rn*eslope/(eslope+psycho)*1e6 / rlat_ht * 100/rho_w)*10
  return(et_max)
}

4. calc_daily_watbal Function (Daily water-balance Function)

Description

Computes the daily water-balance update using the ratio of actual ET (Ea) to potential ET (Ep), given a current water availability and total soil capacities.

Input variable	Output variable	Data source
- soilcp: Field capacity water-holding capacity of the soil (in mm) - soilsat: Additional soil water capacity from field capacity up to saturation (in mm) - avail: Current available soil water content (in mm) - rain: Daily rainfall (in mm/day) - pet: Daily potential evapotranspiration (in mm/day) - cropfc: Crop factor, a multiplier adjusting evapotranspiration for specific crop conditions (default is 1)	- AVAIL: The updated soil water availability (mm) after accounting for ET, rain, and water movement - DEMAND: The evapotranspiration demand (mm), i.e., the amount of water evaporated and transpired - ERATIO: The ratio of actual ET to potential ET (unitless), showing how much of the atmospheric demand was realized given water availability - LOGGING: Water present above field capacity but below saturation (mm), representing temporarily stored water beyond optimal capacity - RUNOFF: Excess water beyond saturation capacity (mm) that leaves the soil system, typically as surface runoff.	- CHIRPS - soilcap_calc - calc_pet

code

# Author: Peter Steward, Fabio Castro, Julian Rameriz
# Date: 2025-03-15
calc_daily_watbal <- function(
    soilcp,
    soilsat,
    avail,
    rain,
    pet,
    cropfc=1
){
  # e/a to e/p subfunction
  eabyep_subfun <- function(soilcp, avail){
    if (soilcp == 0) return(0)
    percwt <- avail/soilcp * 100
    percwt <- max(1, min(100, percwt))
    eratio <- min(percwt/(97 - 3.868*sqrt(soilcp)), 1)
    return(eratio)
  }
  
  # Clamp to field capacity if needed
  new_avail <- min(soilcp, avail)
  
  # Calculate fraction of PET that can be supplied
  eratio <- eabyep_subfun(soilcp, new_avail)
  demand <- eratio*cropfc*pet
  
  # After water usage + new rainfall
  result  <- new_avail + rain - demand
  
  # Water above field capacity
  logging <- result - soilcp
  logging <- ifelse(logging < 0, 0, logging)
  logging <- ifelse(logging > soilsat, soilsat, logging)
  
  # New available water
  new_avail <- min(soilcp, result)
  new_avail <- ifelse(new_avail < 0, 0, new_avail)
  
  # Runoff is beyond saturation
  runoff <- result - (soilcp + logging)
  runoff <- ifelse(runoff < 0, 0, runoff)
  
  data.frame(
    AVAIL   = new_avail,
    DEMAND  = demand,
    ERATIO  = eratio,
    LOGGING = logging,
    RUNOFF  = runoff
  )
}

eabyep_subfun Subfunction

Aspect	Description
Inputs	soilcp and water available capped at field capacity.
Output	Ratio of actual to potential evapotranspiration (eratio), value between 0 and 1.
Function	Empirical model reducing ET proportionally to soil moisture stress.

---