phydro: le is not equal to aet

[fabern]

Using Phydro, the water holding capacity is incorrectly taking into account for the output le.
This is linked to the correction happening when the bucket is empty:

rsofun/src/waterbal_splash.mod.f90

Lines 119 to 129 in 21f4e98

	else if (tile(lu)%soil%phy%wcont < 0.0) then
	! -----------------------------------
	! Bucket is empty
	! -----------------------------------
	! set soil moisture to zero
	tile_fluxes(lu)%canopy%daet = tile_fluxes(lu)%canopy%daet + tile(lu)%soil%phy%wcont
	tile(lu)%soil%phy%wcont = 0.0
	tile_fluxes(lu)%canopy%dro = 0.0
	tile_fluxes(lu)%canopy%dfleach = 0.0
	else

In the above Fortran code, actual transpiration in units of mm (aet = daet) are corrected, but not the actual transpiration in units of energy (le = daet_e).

With Francesco we suggest to add another line that does the same correction for daet_e, too.

• TODO: What still needs to be considered is, if the correction should be applied to the total evapotranspiration (daet) or also the its two constituents: daet_soil and daet_canopy.

[fabern]

Minimum working example (MWE) follows. Below code runs the simulation for FR-Pue for the PHYDRO model for three different values of whc. It should result in different le, corresponding to the differences in aet.

library(rsofun) # 
sofun_4.4.1
library(dplyr)
library(ggplot2)

# define model parameter values from previous
params_modl <- list(
  kphio              = 0.04998,    # setup ORG in Stocker et al. 2020 GMD
  kphio_par_a        = 0.0,        # set to zero to disable temperature-dependence of kphio
  kphio_par_b        = 1.0,
  soilm_thetastar    = 0.6 * 240,  # to recover old setup with soil moisture stress
  soilm_betao        = 0.0,
  beta_unitcostratio = 146.0,
  rd_to_vcmax        = 0.014,      # value from Atkin et al. 2015 for C3 herbaceous
  tau_acclim         = 30.0,
  kc_jmax            = 0.41#,
  # whc                = driver$site_info[[1]]$whc
)


#drivers_for_debugging_PHYDRO <- p_model_drivers # NOT YET UPDATED FOR PHDRO
drivers_for_debugging_PHYDRO <- readRDS(file = here::here("data/p_model_drivers_newformat.rds"))
drivers_for_debugging_PHYDRO$params_siml[[1]]$use_gs     <- TRUE

# Run 3 simulations with different WHC:
output_use_gs <- bind_rows(
  rsofun::runread_pmodel_f(
    drivers_for_debugging_PHYDRO |> tidyr::unnest(site_info) |> mutate(whc = 432) |> tidyr::nest(site_info = !c(sitename, params_siml, starts_with("forcing"))),
    par = purrr::assign_in(params_modl, "whc", 432)
  ) |> mutate(sitename = paste0(sitename, "_432mm")),
  rsofun::runread_pmodel_f(
    drivers_for_debugging_PHYDRO |> tidyr::unnest(site_info) |> mutate(whc = 5) |> tidyr::nest(site_info = !c(sitename, params_siml, starts_with("forcing"))),
    par = purrr::assign_in(params_modl, "whc", 5)
  ) |> mutate(sitename = paste0(sitename, "_5mm")),
  rsofun::runread_pmodel_f(
    drivers_for_debugging_PHYDRO |> tidyr::unnest(site_info) |> mutate(whc = 5000) |> tidyr::nest(site_info = !c(sitename, params_siml, starts_with("forcing"))),
    par = purrr::assign_in(params_modl, "whc", 5000)
  ) |> mutate(sitename = paste0(sitename, "_5000mm"))
)

# Plot:
output_use_gs |> 
  tidyr::unnest(data) |> select(-site_info) |>
  filter(date < "2012-01-01") |>
  select(sitename, date, gpp, aet, le, pet) |>
  tidyr::pivot_longer(!c(sitename,date)) %>%
  ggplot(data = ., mapping=aes(x=date, y=value, color=sitename, linetype=sitename)) +
  geom_line() +
  facet_grid(name~., scales = "free_y") +
  theme_bw()

[fabern]

A test case for the test suite can be to check that aet and le correspond to each other:

# Test for the test suite:
# remotes::install_github("geco-bern/cwd")
# library(cwd)
# convert_et(le, temp, pressure)
convert_et <- function (et_e, tc, patm, return_df = FALSE) {
  # TODO:
  return(et_e)
}
output_use_gs |>
  tidyr::unnest(data) |>
  select(sitename, date, aet, le) |>
  group_by(sitename) |>
  summarise(#identical(aet, aet),
            identical(aet, convert_et(le))) # with the correct function, they should be the same
# IN ABOVE TEST WE NEED TO REQUIRE TRUE

[stineb]

I think the root problem is somewhere else. soilmstress is never below 1.0. I'm looking at it now.

[stineb]

Should be resolved now with 21f4e98.

[fabern]

* [ ]  TODO: What still needs to be considered is, if the correction should be applied to the total evapotranspiration (`daet`) or also the its two constituents: `daet_soil` and `daet_canopy`.

See PR #221 : it makes a suggestion at a correction. However, it is still unsatisfying since daet_canopy can get negative (FR-Pue with WHC of 5mm, see test-model-runs.R)

[fabern]

This is solved by PR #221.

When AET is reduced due to soil moisture stress, we reduce AET_soil and AET_canopy proportianlly, i.e. maintain the fraction of actual canopy transpiration to total actual (AET_canopy / AET), thereby asserting that AET == AET_soil + AET_canopy.