Coupled vegetation-soil-water model

This is the revised version of the full model script used in class, modified to include a sample VPD calculation directly in the driver data.

Model Overview

This toy soil water balance model consists of the following components:

1. Net radiation calculated from incoming shortwave radiation and estimated longwave radiation
   
2. Bucket soil moisture model to simulate water storage and depletion
   
3. Priestley–Taylor evapotranspiration (ET) to estimate atmospheric water demand
   
4. Runoff generation when soil moisture exceeds field capacity
   
5. Soil water potential derived from soil moisture using a simple inverse relationship (1/x curve)

Installing the phydro model

We will be using the develop branch of the Phydro model here:

https://github.com/jaideep777/phydro/tree/develop

Windows users:

1. Install Rtools

2. Install the devtools package

3. run the following command:

devtools::install_github("jaideep777/phydro", ref="develop")

If still facing issues, you may directly install Windows binaries by following the instructions here:

https://github.com/jaideep777/phydro/issues/12

Linux / Mac users:

You may skip step 1 and follow steps 2 and 3 above.

# devtools::install_github("jaideep777/phydro", ref="develop")
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(rphydro)

1. Environmental Drivers

Milestone 1: Example input data

days <- 1:(365*5)

set.seed(42)

monsoon = rep(0,365)
monsoon[150:270] = 1

rain <- monsoon*0.6*rgamma(length(days), shape = 1.2, scale = 5) + 
        0.06*rgamma(length(days), shape = 1.2, scale = 5)  # mm/day rainfall
shortwave <- 700 - 200*cos(days*2*pi/365)            # W/m2 incoming shortwave
Tair <- -5*cos(days*2*pi/365) + 25 + runif(length(days), -5, 5)                   # air temperature (°C)
esat = Tair |> sapply(rphydro::calc_esat, rphydro::calc_patm(0))
rh = 0.6 + scale(rain, center=F)/2. #relative humidity
rh[rh>0.98] = 0.98
vpd = esat*(1-rh)*1e-3 # VPD in kPa. #vapor pressure deficit
data <- data.frame(days, rain, shortwave, Tair, esat, vpd, rh)

data |> 
  pivot_longer(-days) |> 
  ggplot(aes(x=days, y=value)) +
  geom_line() + 
  facet_wrap(~name, scales="free")

2. Model parameters

The model represents soil as a single-layer storage reservoir defined by two moisture thresholds.

• Field capacity (θ_fc): maximum water the soil can retain after drainage
  
• Wilting point (θ_wp): minimum soil moisture at which plants can extract water

When soil moisture approaches the wilting point, plants experience water stress and evapotranspiration is limited.

# Physical constants
sigma <- 5.67e-8   # Stefan-Boltzmann constant (W m-2 K-4)
albedo <- 0.23     # surface albedo
alpha_PT <- 1.26   # Priestley-Taylor coefficient
lambda <- 2.257e6   # latent heat of vaporization (J/kg)
gamma <- 0.066     # psychrometric constant (kPa/C)

# ---------------------------
# Soil parameters
# ---------------------------

soil_depth <- 150              # mm storage capacity
theta_fc <- soil_depth
theta_wp <- 100
theta_init <- 120              # initial soil water
k_psi <- -2                 # scaling constant for potential

3. Supporting funcitions

Milestone 2: Vapour presure and emmissivity functions

Saturation vapour pressure (SVP) increases nonlinearly with temperature following the Clausius–Clapeyron relationship.

Atmospheric emissivity determines what fraction of the surface’s outgoing longwave radiation is re-emitted back downward by the atmosphere. The Brutsaert approximation used here relates emissivity to the actual vapour pressure near the surface — moister air absorbs and re-radiates more longwave radiation, reducing net radiative cooling of the surface.

# Saturation vapour pressure vs celsius temp (kPa)
svp = function(Tc){
  es <- 0.6108 * exp((17.27*Tc)/(Tc+237.3))
  es
}

# saturation vapor pressure slope (kPa/K)
delta_svp <- function(T){
  es <- svp(T)
  delta <- (4098 * es)/((T+237.3)^2)
  return(delta)
}

# estimate atmospheric emissivity (Brutsaert approximation: https://www.nature.com/articles/s41598-023-40499-6#Sec7) 
# ea is actual vapour pressure in kPa
# paper says ea is in hPa but that gives unrealistic values, likely should be kPa
atm_emissivity <- function(T, ea){
  T_k <- T + 273.15
  ea <- pmax(ea,0)
  eps <- 1.24 * (ea/T_k)^(1/7)  # requires ea in hPa
  pmin(eps,1)
}


print(
tibble(temp = 1:45) |> 
  mutate(vpd = 0.6) |> # vpd in kPa
  mutate(saturation_vp = svp(temp)) |> 
  mutate(delta_svpvst = delta_svp(temp)) |> 
  mutate(emissivity = atm_emissivity(temp, saturation_vp - vpd)) |> 
  pivot_longer(-temp) |> 
  ggplot(aes(x=temp, y=value)) + 
  geom_line() + 
  facet_wrap(~name, scales = "free")
)

Milestone 3: Net radiation calculation

Net radiation (Rn) is the primary energy driver of the system. It is calculated as the balance between absorbed shortwave radiation (incoming solar minus reflected) and net longwave radiation (surface emission minus atmospheric back-radiation). Net radiation is positive during the day when solar input dominates, and can be negative at night or under low-humidity conditions when the surface loses more longwave than it receives back.

# net radiation calculation
# shortwave in W/m2, vpd in kPa, Tc in degC
net_radiation <- function(shortwave, vpd, Tc){

  es <- 0.6108 * exp((17.27*Tc)/(Tc+237.3))
  ea = es - vpd

  # Net shortwave
  Rns <- (1 - albedo) * shortwave

  # Longwave radiation
  T_k <- Tc + 273.15
  eps_a <- atm_emissivity(Tc, ea)

  Rl_down <- eps_a * sigma * T_k^4
  Rl_up <- sigma * T_k^4

  Rnl <- Rl_up - Rl_down

  Rn <- Rns - Rnl
  return(Rn)
}

Milestone 4: Implement Priestly-Taylor ET

Evapotranspiration (Priestley–Taylor)

Evapotranspiration (ET) describes the combined water loss from soil evaporation and plant transpiration. Potential ET is estimated using the Priestley–Taylor approach. In the model, actual ET (AET) is further constrained by soil moisture — when soil water approaches the wilting point, AET is reduced proportionally below the potential rate.

phydro_ET <- function(fapar, Rn, tc, vpd_kpa, psi_soil){
  kphio = 0.04;        # quantum yield efficiency
  ppfd = max(Rn*2, 0);          # umol/m2/s
  co2_ppm  = 400;          # ppm
  elv  = 0;            # m.a.s.l.
  rdark = 0.015;
  vwind = 3;
  netrad = ppfd/2
  pa  = calc_patm(elv);            # m.a.s.l.

  vpd_pa = vpd_kpa*1000

  par_cost = list(alpha=0.1, gamma=1);
  par_plant = list(conductivity=3e-17, psi50=-2, b=2);
  options = list(gs_method = "GS_IGF", 
                et_method = "ET_DIFFUSION",
                ftemp_vj_method = "FV_kumarathunge19",
                ftemp_rd_method = "FR_heskel16",
                ftemp_br_method = "FB_atkin15",
                scale_alpha = F)  
  
  res = rphydro_analytical(tc, tc, ppfd, netrad, vpd_pa, co2_ppm, pa, fapar, kphio, psi_soil, rdark, vwind, par_plant, par_cost, options)

  list(
    et = 1.6*res$gs*(vpd_pa/pa)*18.015*1e-3*86400,  # mm / day
    a = res$a*1e-6*12.01*86400   # gC / m2 / day
  )
}

# Priestley–Taylor ET
priestley_taylor_ET <- function(Rn, Tc){

  delta <- delta_svp(Tc) # (kPa/K)

  ET <- alpha_PT * (delta/(delta + gamma)) * Rn

  # convert W/m2 -> mm/day
  ET_mm <- ET * 86400 / lambda * 1e-3 * 1000  # J/s/m2 * s/day * kg/J * m3/kg * mm/m

  return(max(0, ET_mm))
}

# Soil moisture -> soil water potential
soil_water_potential <- function(theta){

  psi <- k_psi / theta   # simple 1/x curve

  return(psi)
}

4. Model Run

Milestone 5: Implement water bucket model

The bucket model treats the soil column as a single well-mixed reservoir. At each daily timestep, water enters via rainfall and leaves via actual evapotranspiration (AET) and runoff/drainage.

The LAI (Leaf Area Index) evolves dynamically as a simple relaxation toward a target value that depends on available energy (Rn), representing the tendency of vegetation to grow more leaf area under favourable conditions and shed it otherwise.

The phydro ET model replaces the simpler Priestley–Taylor formulation by explicitly modelling plant hydraulic limitation — as soil water potential becomes more negative, stomatal conductance is reduced, limiting both transpiration and photosynthesis (GPP). This couples the carbon and water cycles within the same framework.

run_and_plot = function(){
  n <- nrow(data)
  
  soil_water <- numeric(n)
  PET <- numeric(n)
  AET <- numeric(n)
  runoff <- numeric(n)
  Rn <- numeric(n)
  psi <- numeric(n)
  gpp <- numeric(n)
  lai <- numeric(n)
  
  soil_water[1] <- theta_init
  psi[1] <- soil_water_potential(theta_init)
  lai[1] <- 1
  
  # ---------------------------
  # Simulation loop
  # ---------------------------
  for(i in 2:n){
    # Net radiation
    Rn[i] <- net_radiation(data$shortwave[i], data$vpd[i], data$Tair[i])
    
    # LAI model
    lai[i] <- lai[i-1]*0.95 + 0.05*(1 + (Rn[i] > 400)*4)
    
    # # Potential ET (Priestley–Taylor)
    # PET[i] <- priestley_taylor_ET(Rn[i], data$Tair[i])
    
    # Add rainfall
    soil_water[i] <- soil_water[i-1] + data$rain[i]
    
    # # -------------------------
    # # Soil moisture stress
    # # -------------------------
    
    # stress <- (soil_water[i] - theta_wp) / (theta_fc - theta_wp)
    # stress <- max(0, min(1, stress))   # constrain 0–1
    
    # AET[i] <- PET[i] * stress
    # # -------------------------
    # # Direct AET calculation from phydro
    # # -------------------------
    res <- phydro_ET(fapar=1-exp(-0.5*lai[i]), Rn[i], data$Tair[i], data$vpd[i], psi[i-1])
    AET[i] <- res$et
    gpp[i] <- res$a
    
    # Do not extract water below wilting point
    available_water <- soil_water[i] - theta_wp
    AET[i] <- min(AET[i], available_water)
    
    # Update soil storage
    soil_water[i] <- soil_water[i] - AET[i]
    
    # -------------------------
    # Runoff / drainage
    # -------------------------
    if (soil_water[i] > theta_fc){
      
      # optional gravitational drainage
      runoff[i] <- soil_water[i] - theta_fc
      soil_water[i] <- theta_fc
      
    } else {
      
      runoff[i] <- 0
      
    }
    
    # -------------------------
    # Soil water potential
    # -------------------------
    
    psi[i] <- soil_water_potential(soil_water[i]-theta_wp)
    
  }
  
  # ---------------------------
  # Combine results
  # ---------------------------
  
  results <- data.frame(
    day = data$days,
    rainfall = data$rain,
    # shortwave = data$shortwave,
    lai = lai,
    temperature = data$Tair,
    net_radiation = Rn,
    # pet = PET,
    aet = AET,
    gpp = gpp,
    soil_water = soil_water,
    runoff = runoff,
    soil_water_potential = psi
  ) |> slice(-1)
  
  
  # ---------------------------
  # Plot results
  # ---------------------------
  print(
    results |> 
      mutate(day= as.Date("2021-01-01")+day) |> 
      # select(day, soil_water, evapotranspiration, soil_water_potential) |> 
      pivot_longer(-day) |> 
      ggplot(aes(x=day, y=value)) +
      geom_line() +
      facet_wrap(~name, scales="free")
  )
  
  results |> colMeans() |> print()
}

#default run
run_and_plot()

##                  day             rainfall                  lai 
##          913.5000000            1.4172720            2.7716702 
##          temperature        net_radiation                  aet 
##           25.0662673          369.8263198            0.7495444 
##                  gpp           soil_water               runoff 
##            4.0617421          120.8889709            0.6748786 
## soil_water_potential 
##           -0.8316063

Appendix: Tutorial milestones for submission

1. Change the environmental drivers to represent moist conditions with very low seasonality and run the model. Comment on the output and how it differs from the demo output shown above.

days <- 1:(365*5)

set.seed(42)

monsoon = rep(0,365)
monsoon[1:365] = 1

rain <- monsoon*0.2*rgamma(length(days), shape = 1.2, scale = 5) + 
  0.06*rgamma(length(days), shape = 1.2, scale = 5)  # mm/day rainfall
shortwave <- 700 - 200*cos(days*2*pi/365)            # W/m2 incoming shortwave
Tair <- -5*cos(days*2*pi/365) + 25 + runif(length(days), -5, 5)                   # air temperature (°C)
esat = Tair |> sapply(rphydro::calc_esat, rphydro::calc_patm(0))
rh = 0.6 + scale(rain, center=F)/2
rh[rh>0.98] = 0.98
vpd = esat*(1-rh)*1e-3 # VPD in kPa
data <- data.frame(days, rain, shortwave, Tair, esat, vpd, rh)

data |> 
  pivot_longer(-days) |> 
  ggplot(aes(x=days, y=value)) +
  geom_line() + 
  facet_wrap(~name, scales="free")

run_and_plot()

##                  day             rainfall                  lai 
##          913.5000000            1.5396407            2.8133616 
##          temperature        net_radiation                  aet 
##           25.0662673          376.4329665            0.8347509 
##                  gpp           soil_water               runoff 
##            4.8967941          148.4710701            0.6884424 
## soil_water_potential 
##           -0.0415223

Hint: Removing the monsoon signal and distributing rainfall more uniformly throughout the year reduces the amplitude of soil moisture fluctuations. With soil moisture remaining closer to field capacity year-round, water stress is rarely triggered, so AET tracks PET more closely and GPP remains relatively stable across seasons. The absence of a pronounced dry season removes the main driver of stomatal closure and plant stress.

2. Change the plant parameters (plant_par object in the phydro_ET function) to represent a drought-avoiding plant. Such a plant will have a less negative psi50, say -0.5 or -0.75. Comment on the outcomes in comparison with the default scenario.

#{Restore environmental drivers to original values}
days <- 1:(365*5)

set.seed(42)

monsoon = rep(0,365)
monsoon[150:270] = 1

rain <- monsoon*0.6*rgamma(length(days), shape = 1.2, scale = 5) + 
        0.06*rgamma(length(days), shape = 1.2, scale = 5)  # mm/day rainfall
shortwave <- 700 - 200*cos(days*2*pi/365)            # W/m2 incoming shortwave
Tair <- -5*cos(days*2*pi/365) + 25 + runif(length(days), -5, 5)                   # air temperature (°C)
esat = Tair |> sapply(rphydro::calc_esat, rphydro::calc_patm(0))
rh = 0.6 + scale(rain, center=F)/2
rh[rh>0.98] = 0.98
vpd = esat*(1-rh)*1e-3 # VPD in kPa
data <- data.frame(days, rain, shortwave, Tair, esat, vpd, rh)

data |> 
  pivot_longer(-days) |> 
  ggplot(aes(x=days, y=value)) +
  geom_line() + 
  facet_wrap(~name, scales="free")

#change the plant p50 to -0.50
par_plant = list(conductivity=3e-17, psi50=-0.5, b=2);

run_and_plot()

##                  day             rainfall                  lai 
##          913.5000000            1.4172720            2.7716702 
##          temperature        net_radiation                  aet 
##           25.0662673          369.8263198            0.7495444 
##                  gpp           soil_water               runoff 
##            4.0617421          120.8889709            0.6748786 
## soil_water_potential 
##           -0.8316063

Hint: A less negative psi50 means stomata begin closing at a much higher (less negative) soil water potential, thus, the plant avoids hydraulic failure by shutting down early. Compared to the default (psi50 = −2), this plant will reduce transpiration and GPP sooner as the soil dries, maintaining higher soil moisture but at the cost of reduced carbon assimilation during dry periods.

3. Increase the water storage capacity of the soil by increasing the soil depth and run the simulation. Comment on the outcomes and compare them with the default run.

#Restore plant paramters and set soil dept
par_plant = list(conductivity=3e-17, psi50=-2, b=2);

soil_depth <- 300              # mm storage capacity
theta_fc <- soil_depth

run_and_plot()

##                  day             rainfall                  lai 
##          913.5000000            1.4172720            2.7716702 
##          temperature        net_radiation                  aet 
##           25.0662673          369.8263198            0.9388605 
##                  gpp           soil_water               runoff 
##            4.4749876          238.8297817            0.4036479 
## soil_water_potential 
##           -0.1525430

Hint: Doubling soil depth increases the total water reservoir available to plants. This acts as a buffer against drought as the same rainfall events recharge a proportionally smaller fraction of capacity, and the soil sustains plant water supply for longer into dry periods. The result is reduced inter-annual variability in AET and GPP, and less frequent triggering of the wilting-point stress threshold.