Faulty pattern of water absorption with deeper roots

[MichaKoe]

I am modelling the water budget in 2021-2022 of a forest stand in southern Hesse. 2021 was rather "wet" while 2022 was a very dry year.

I have a 4m deep bottom profile and set the parameters:

maxrootdepth <- -4
betaroot<-0.99
root_method<- "linear"
psiini <- (-15) 

The rooting patterns extracted from model input looks like this:

The resulting water budget and matric potentials:

It looks like the model allocates the roots mainly below 3m depth and practically hardly absorbs any water from the topsoil. Probably there is a sorting error in the transfer of the rooting patterns from preprocessing to the model.

Edit: it does not change the picture if I am setting bypar=0 Edit2: I am working with 15 min precipitation interval

[MichaKoe]

Keeping all parameters the same but changing maxrootdepth to -2 gives the following results: 2m

[MichaKoe]

I made a reproducible exampe with data from within the package only varying maximum rooting depth:

library(LWFBrook90R)
library(data.table)
library(ggplot2)

par(mfrow=c(2,2))
data("slb1_meteo")

xx<- -4

out<-lapply(c(-1,-2,-3,-4), function(xx){

  options_b90 <- set_optionsLWFB90()
  param_b90 <- set_paramLWFB90()
  param_b90$maxrootdepth<- xx
  param_b90$betaroot<-0.99
  soil <- cbind(slb1_soil, hydpar_wessolek_tab(texture = slb1_soil$texture))
  options_b90$startdate
  options_b90$enddate<-as.Date("2002-12-31")

  # create some extra soil below
  extrasoil<-soil[rep(nrow(soil), 30),]
  soil<-soil[-nrow(soil),]
  extrasoil$upper<-seq(min(soil$lower),-3.9,length.out=nrow(extrasoil))
  extrasoil$horizon<-"fake"
  soil<-rbind(soil[-(nrow(soil)),],extrasoil)
  soil$lower<-c(soil$upper[-1],-4)

  b90res <- run_LWFB90(options_b90 = options_b90,
                       param_b90 = param_b90,
                       climate = slb1_meteo,
                       soil = soil)

  output <- set_outputLWFB90()
  output[,] <- 0 # reset output
  output["Swat", "Day"] <- 1
  output[c("Evap", "Flow","Abov"), "Mon"] <- 1

  results<-process_outputs_LWFB90(b90res,
                                  selection = output,
                                  prec_interval = options_b90$prec_interval)

  soil_out<-b90res$model_input$param_b90$soil_nodes
  soil_out$nl<-soil_out$layer

  plot(y=soil_out$midpoint,soil_out$rootden,type="b",
       main=paste0("rooting depth",xx))

  mpot_mean <- merge(results$SWATDAY.ASC, soil_out[,list(nl,upper,lower)],
                     by = c("nl"))

  mpot_mean[,date:=as.Date(paste0(yr,"-",mo,"-",da))]

  # interpolate from discrete soil depths to "continuous" cm steps by day of year
  mpot_mean <- mpot_mean[, approx(lower * (-100), psimi,
                                  xout = (0:400),
                                  rule = 2), by = list( doy)]

  breaks_psimi <- -c(2000,1000,250,225,200,175,150,125,100,75,50,25,10,5,0)
  mpot_mean[,y_discr := cut(y, breaks = breaks_psimi,)]

  colfunc <- colorRampPalette(rev(c("blue", "white","red")))

  p_mpot <- ggplot(mpot_mean, aes(doy, x)) +
    geom_raster(aes(fill = y_discr)) +
    geom_hline(aes(yintercept = param_b90$maxrootdepth*-100),
               linetype = "dashed", color = "white",
               linewidth = 0.4)+
    scale_y_reverse(sec.axis = sec_axis(trans = ~.*1), expand = c(0,0)) +
    scale_x_continuous(expand = c(0,0)) +
    ggtitle(paste0("routing depth: ",xx ))+
    scale_fill_manual(values = colfunc(length(unique(mpot_mean$y_discr))))+
    labs(x = "Day of year", y = "Soil depth [cm]", fill =
           expression(paste(Psi, " [kPa]")))

    return(p_mpot)
})
names(out)<-paste0("rd_",c(-1,-2,-3,-4))

print(gridExtra::grid.arrange(out$`rd_-1`,
                              out$`rd_-2`,
                              out$`rd_-3`,
                              out$`rd_-4`, nrow = 2))

This gives the following output:

[pschmidtwalter]

I found the cause. The root density depth distribution is transferred correctly to the Fortran subroutine. However, the root-growth subroutine then calculates a daily growing root distribution, depending on the age of the trees, approaching the provided final root distribution. An unfortunate combination of tree age, rgrorate, and rgroper and maxrootdepth caused roots growing into the lowest layers. In these lowest layers, where the root density should be 0 at the current age (100), something bad happens.

A quick work-around that prevents root growth is to set param_b90$rgroper = 0 .

This will also be the new default value for this parameter.

[MichaKoe]

I took a closer look and to be honest, I don't quite understand the root growth part of the model. Hammel and Kennel say: Depth growth is linear during some “juvenile phase”, going from a starting depth to a maximum depth. The slope of this linear function is thus the growth rate in m/yr. Now the standard parameters that LWFBrook90R passes to the underlying LWFBrook90 are the following:

• “Initial root depth”initrdep=0.25 m (Note: I think this default value is wrongly positive in LWFBrook90Rs standard paramters as root depths are actually always negative in Brook90)

• the maximum root depth maxrootdepth= -1.5 m (for Hammel and Kennel this is given by the input root distribution, for LWFBrook90R this is used to generate the input root distribution via "make_rootden" )

• “Period of net root growth” rgroper = 30 years (juvenile phase)

• the "Vertical root growth rate" rgrorate = 0.03 m per year

This last parameter is unclear to me, because if the vertical rate should be meant by this, then it would exist twice in the model: it already results from (maximum root depth - initial root depth) / growth duration.

Example: if you assume this standard vertical growth rate, then at the end of the active root growth for a 30 year old stand you get 90 cm depth + 25 cm initial root depth = 1.15 m maximum root depth and not the required -1.5 m maxrootdepth of the input root distribution.

Hammel and Kennel also intended to model lateral root growth. According to them the roots grow laterally, i.e. increase their “relative root length density” during growth in every next layer reached. This value is calculated based on the total root length and the relative root distribution. But since they also provide initial (initrlen) and final values (maxrlen) also no growth rate (in m/m2/yr or m/yr) is needed here. They even say the “rate parameter” for lateral growth in each depth is “individually adjusted” such that the given final root distribution is reached until the end of the period of root growth. So my question is: Why is there a need for a growth rate (vertical or lateral)?

---

Further: I modeled without root growth. My stand was 100 years old, I even set

  param_b90$initrdep<- xx
  param_b90$maxrootdepth<- xx
  param_b90$initrlen <- param_b90$maxrlen 

Therefore the mentioned problem "roots never growing in to the lowest layers" should no apply to my cases because I tell the model that they are already there!

To conclude: I am not getting it.

[pschmidtwalter]

A quick work-around that prevents root growth is to set param_b90$rgroper = 0 .

This will also be the new default value for this parameter.

As of version 0.5.2 the input-parameter param_b90$rgroper is set to 0 by default, preventing root growth and forcing the model to use the final root length density depth distribution in the calculations. So, this needs to be actively switched on by the user if age dependend root growth should be used, e.g. for juvenile forest stands. However, layer-wise water uptake (layer_output$tran) should be carefully investigated when using root growth. I will try to fix the problems associated with the rootgrowth-module in a later version.