In emf-creaf/medfateland: Mediterranean Landscape Simulation

knitr::opts_chunk$set(echo = TRUE)
library(medfate)
library(medfateland)
library(meteoland)
library(ggplot2)

Introduction

The aim of this vignette is to illustrate the sensitivity of watershed ecohydrological simulations (e.g. via spwb_land()) to variations in three watershed parameters for sub-model TETIS. These parameters are scaling factors of hydraulic conductivities for soil vertical fluxes, subsurface lateral fluxes and groundwater lateral fluxes, respectively. They are referred to as R_localflow, R_interflow and R_baseflow, respectively in the package.

Preparation

Here we load a small example watershed included with the package, that can be used to understand the inputs required:

data("example_watershed_burnin")
example_watershed_burnin

The watershed state is the result of a four-year burn-in period, in order to have aquifer and soil moisture at a (dynamic) equilibrium. Another spatial input is needed to describe the grid topology, which in our case is an object of class SpatRaster from package terra:

r <-terra::rast(xmin = 401380, ymin = 4671820, xmax = 402880, ymax = 4672620, 
                nrow = 8, ncol = 15, crs = "epsg:32631")
r

Finally, we load an example weather data set to be used as climate forcing, as well as the usual species parameter table:

# Load example weather dataset
data("examplemeteo")
# Set simulation period
dates <- seq(as.Date("2001-01-01"), as.Date("2001-03-31"), by="day")
# Load default medfate parameters
data("SpParamsMED")

Initial values for scaling factors

Watershed simulations have overall control parameters. Notably, the user needs to decide which sub-model will be used for lateral water transfer processes, in this case:

ws_control <- default_watershed_control("tetis")

The default parameterization for the three scaling factors is:

ws_control$tetis_parameters

Simulations

Since watershed simulations are time consuming, here we will only perform four simulations, i.e. a baseline plus a modification of each of the three scaling factors, one at a time.

Baseline configuration For the baseline configuration, we set the three scaling factors to one:

ws_control$tetis_parameters$R_localflow <- 1
ws_control$tetis_parameters$R_interflow <- 1
ws_control$tetis_parameters$R_baseflow <- 1

Now we launch the simulation. Focus on the values (in mm) predicted for aquifer exfiltration, saturation excess and, finally, the sum of the two sources as export runoff.

res_0 <- spwb_land(r, example_watershed_burnin, SpParamsMED, examplemeteo, 
                   dates = dates, summary_frequency = "month", summary_blocks = "WaterBalance",
                   watershed_control = ws_control, progress = FALSE)

In the following we compare these values to those issued from simulations with altered scaling factors.

Increasing local vertical conductivity

We begin our sensitivity analysis by increasing the vertical conductivity five times:

ws_control$tetis_parameters$R_localflow <- 5
ws_control$tetis_parameters$R_interflow <- 1
ws_control$tetis_parameters$R_baseflow <- 1
res_L5 <- spwb_land(r, example_watershed_burnin, SpParamsMED, examplemeteo, 
                     dates = dates, summary_frequency = "month",summary_blocks = "WaterBalance",
                     watershed_control = ws_control, progress = FALSE)

Increasing R_localflow has the effect of increasing the drainage from soils to the aquifer (hence, decreasing soil moisture). It also increases capillarity rise due to the higher conductivity. Aquifer exfiltration increases downslope as a consequence of the larger recharge in upslope cells.

Increasing subsurface lateral conductivity

We now focus on the lateral conductivity of subsurface fluxes (interflow):

ws_control$tetis_parameters$R_localflow <- 1
ws_control$tetis_parameters$R_interflow <- 5
ws_control$tetis_parameters$R_baseflow <- 1
res_I5 <- spwb_land(r, example_watershed_burnin, SpParamsMED, examplemeteo, 
                     dates = dates, summary_frequency = "month",summary_blocks = "WaterBalance",
                     watershed_control = ws_control, progress = FALSE)

Increasing lateral conductivity for subsurface fluxes increases the rate of lateral exchange (i.e. the speed of the kinematic wave), resulting in a much larger saturation excess. It also impacts the amount of aquifer exfiltration.

Increasing groundwater lateral conductivity Let's now perform a simulation with an increased groundwater conductivity:

ws_control$tetis_parameters$R_localflow <- 1
ws_control$tetis_parameters$R_interflow <- 1
ws_control$tetis_parameters$R_baseflow <- 5
res_B5 <- spwb_land(r, example_watershed_burnin, SpParamsMED, examplemeteo, 
                     dates = dates, summary_frequency = "month",summary_blocks = "WaterBalance",
                     watershed_control = ws_control, progress = FALSE)

Increasing the speed of groundwater fluxes leads to a much shallower water table downslope, significantly increasing the amount of exfiltration.

Comparing predicted fluxes

Flux rate series

df1<- res_0$watershed_balance
df1$Scenario <- "Baseline"
df2<- res_L5$watershed_balance
df2$Scenario <- "R_local x 5"
df3<- res_I5$watershed_balance
df3$Scenario <- "R_inter x 5"
df4<- res_B5$watershed_balance
df4$Scenario <- "R_base x 5"
df<- rbind(df1,df2,df3, df4)
df$Scenario <- factor(df$Scenario, levels = c("Baseline", "R_local x 5", "R_inter x 5", "R_base x 5"))

Drainage

ggplot(df)+
  geom_line(aes(x=dates, y =DeepDrainage, col = Scenario, linetype = Scenario))+
  ylab("Deep drainage (mm)")+
  theme_bw()

Aquifer exfiltration

ggplot(df)+
  geom_line(aes(x=dates, y =AquiferExfiltration, col = Scenario, linetype = Scenario))+
  ylab("Aquifer exfiltration (mm)")+
  theme_bw()

Flux rate maps

Drainage to aquifer

g1 <- plot_summary(res_0$sf, "DeepDrainage", date = "2001-01-01", r = r, limits = c(0, 100))+labs(title = "Baseline")
g2 <- plot_summary(res_L5$sf, "DeepDrainage", date = "2001-01-01", r = r, limits = c(0, 100))+labs(title = "R_local x 5")
g3 <- plot_summary(res_I5$sf, "DeepDrainage", date = "2001-01-01", r = r, limits = c(0, 100))+labs(title = "R_inter x 5")
g4 <- plot_summary(res_B5$sf, "DeepDrainage", date = "2001-01-01", r = r, limits = c(0, 100))+labs(title = "R_base x 5")
cowplot::plot_grid(g1, g2, g3, g4 , ncol = 1, nrow = 4)

Sub-surface balance

g1 <- plot_summary(res_0$sf, "InterflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(limits = c(-100, 450), mid="gray", na.value = NA)+
  labs(title = "Baseline")
g2 <- plot_summary(res_L5$sf, "InterflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(limits = c(-100, 450), mid="gray", na.value = NA)+
  labs(title = "R_local x 5")
g3 <- plot_summary(res_I5$sf, "InterflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(limits = c(-100, 450), mid="gray", na.value = NA)+
  labs(title = "R_inter x 5")
g4 <- plot_summary(res_B5$sf, "InterflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(limits = c(-100, 450), mid="gray", na.value = NA)+
  labs(title = "R_base x 5")
cowplot::plot_grid(g1, g2, g3, g4 , ncol = 1, nrow = 4)

Groundwater balance

g1 <- plot_summary(res_0$sf, "BaseflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(mid="gray", limits = c(-50, 250), na.value = NA)+
  labs(title = "Baseline")
g2 <- plot_summary(res_L5$sf, "BaseflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(mid="gray", limits = c(-50, 250), na.value = NA)+
  labs(title = "R_local x 5")
g3 <- plot_summary(res_I5$sf, "BaseflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(mid="gray", limits = c(-50, 250), na.value = NA)+
  labs(title = "R_inter x 5")
g4 <- plot_summary(res_B5$sf, "BaseflowBalance", date = "2001-03-01", r = r)+
  scale_fill_gradient2(mid="gray", limits = c(-50, 250), na.value = NA)+
  labs(title = "R_base x 5")
cowplot::plot_grid(g1, g2, g3, g4 , ncol = 1, nrow = 4)