In emf-creaf/medfateland: Mediterranean Landscape Simulation

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

Aim

The aim of this vignette is to illustrate how to use medfateland (v. r packageVersion("medfateland")) to carry out simulations of forest function and dynamics on a set of forest stands while including lateral water transfer processes. This is done using functions spwb_land(), growth_land() and fordyn_land(); which are counterparts of functions spwb(), growth() and fordyn() in package medfate. We will focus here on function spwb_land(), but the other two functions would be used similarly. The same can be said for functions spwb_land_day() and growth_land_day(), which are counterparts of spwb_day() and growth_day(), respectively.

Preparation

Preparing inputs for watershed simulations can be tedious. Two main inputs need to be assembled, described in the following two sections (see also Preparing inputs II: arbitrary locations).

Input sf objects

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

data("example_watershed")
example_watershed

Some of the columns like forest, soil, elevation, or state, were also present in the example for spatially-uncoupled simulations, so we will not repeat them. The following describes additional columns that are relevant here.

Land cover type

Simulations over watersheds normally include different land cover types. These are described in column land_cover_type:

table(example_watershed$land_cover_type)

Local and landscape processes will behave differently depending on the land cover type.

Aquifer and snowpack

Columns aquifer and snowpack are used as state variables to store the water content in the aquifer and snowpack, respectively.

Crop factors

Since the landscape contains agricultural lands, we need to define crop factors, which will determine transpiration flow as a proportion of potential evapotranspiration:

example_watershed$crop_factor = NA
example_watershed$crop_factor[example_watershed$land_cover_type=="agriculture"] = 0.75

Channel network

In large watersheds, the hydrological behavior of the model may not be appropriate because water routing in the river channel is not considered. If a binary column called channel is included in the input, the model will use it to determine the river network, outlets and the time in days to reach them (see function overland_routing for a static analysis of channel routing).

Grid topology

Note that the sf structure does not imply a grid per se. Point geometry is used to describe the central coordinates of grid cells, but does not describe the grid. This means that 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

The r object must have the same coordinate reference system as the sf object. Moreover, each grid cell can contain up to one point of the sf (typically at the cell center). Some grid cells may be empty, though, so that the actual simulations may be done on an incomplete grid. Note that the raster does not contain data, only the topology is needed (to define neighbors and cell sizes, for example). All relevant attribute data is already included in the sf object.

Combining the r and sf objects allows drawing rasterized maps:

plot_variable(example_watershed, variable = "elevation", r = r)

Watershed control options

Analogously to local-scale simulations with medfate, watershed simulations have overall control parameters. Notably, the user needs to decide which sub-model will be used for lateral water transfer processes (a decision similar to choosing the plant transpiration sub-model in medfate), by default "tetis":

ws_control <- default_watershed_control("tetis")

Initialization

Simulation model inputs need to be created for the target watershed before launching simulations. This may be done automatically, though, when calling watershed simulation functions, but in many occasions it is practical to perform this step separately. If we plan to use function spwb_land(), watershed initialization would be as follows:

example_init <- initialize_landscape(example_watershed, SpParams = SpParamsMED,
                                     local_control = defaultControl(soilDomains = "buckets"))
example_init

Here we use function defaultControl() to specify the control parameters for local processes. Function initialize_landscape() makes internal calls to spwbInput() of medfate and defines a column state with the initialized inputs.

At this point is important to learn one option that may speed up calculations. Initialization may be done while simplifying forest structure to the dominant species (see function forest_reduceToDominant() in package medfate). Hence, we can initialize using simplify = TRUE:

example_simplified <- initialize_landscape(example_watershed, SpParams = SpParamsMED,
                                           local_control = defaultControl(),
                                           simplify = TRUE)
example_simplified

For computational reasons, we will keep with this simplified initialization in the next sections.

Carrying out simulations

Launching watershed simulations

To speed up calculations, we call function spwb_land() for a single month:

dates <- seq(as.Date("2001-01-01"), as.Date("2001-01-31"), by="day")
res_ws1 <- spwb_land(r, example_simplified,
                    SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month",
                    watershed_control = ws_control, progress = FALSE)

Although simulations are performed using daily temporal steps, parameter summary_frequency allows storing cell-level results at coarser temporal scales, to reduce the amount of memory in spatial results (see also parameter summary_blocks to decide which outputs are to be kept in summaries).

Structure of simulation outputs

Function spwb_land() and growth_land() return a list with the following elements:

names(res_ws1)

Where sf is an object of class sf, analogous to those of functions *_spatial():

res_ws1$sf

Columns state, aquifer and snowpack contain state variables, whereas summary contains temporal summaries for all cells. Column result is empty in this case, but see below.

The next two elements of the simulation result list, namely watershed_balance and watershed_soil_balance, refer to watershed-level results. For example, watershed_balance contains the daily elements of the water balance at the watershed level, including the amount of water exported in mm in the last column.

head(res_ws1$watershed_balance)

Values of this output data frame are averages across cells in the landscape. Data frame watershed_soil_balance is similar to watershed_balance but focusing on cells that have a soil (i.e. excluding artificial, rock or water land cover). Finally, channel_export_m3s contains the average river flow reaching each channel cell each day and outlet_export_m3s contains the average river flow reaching each outlet cell each day (both in units of cubic meters per second):

head(res_ws1$outlet_export_m3s)

Accessing and plotting cell summaries

Unlike spwb_spatial() where summaries could be arbitrarily generated a posteriori from simulation results, with spwb_land() the summaries are always fixed and embedded with the simulation result. For example, we can inspect the summaries for a given landscape cell using:

res_ws1$sf$summary[[1]]

Additional variables (water balance components, carbon balance components, etc.) can be added to summaries via parameter summary_blocks. Some of these summaries are temporal averages (e.g. state variables), while others are temporal sums (e.g. water or carbon balance components), depending on the variable.

Maps of variable summaries can be drawn from the result of function spwb_land() in a similar way as done for spwb_spatial(). As an example we display a map of the average soil relative water content during the simulated month:

plot_summary(res_ws1$sf, variable = "RWC", date = "2001-01-01", r = r)

Full simulation results for specific cells

The idea of generating summaries arises from the fact that local models can produce a large amount of results, of which only some are of interest at the landscape level. Nevertheless, it is possible to specify those cells for which full daily results are desired. This is done by adding a column result_cell in the input sf object:

# Set request for daily model results in cells number 3 and 9
example_simplified$result_cell <- FALSE
example_simplified$result_cell[c(3,9)] <- TRUE

If we launch the simulations again (omitting progress information):

res_ws1 <- spwb_land(r, example_simplified,
                    SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month",
                    watershed_control = ws_control, progress = FALSE)

We can now retrieve the results of the desired cell, e.g. the third one, in column result of sf:

S <- res_ws1$sf$result[[3]]
class(S)

This object has class spwb and the same structure returned by function spwb() of medfate. Hence, we can inspect daily results using functions shinyplot() or plot(), for example:

plot(S, "SoilRWC")

Continuing a previous simulation

The result of a simulation includes an element state, which stores the state of soil and stand variables at the end of the simulation. This information can be used to perform a new simulation from the point where the first one ended. In order to do so, we need to update the state variables in spatial object with their values at the end of the simulation, using function update_landscape():

example_watershed_mod <- update_landscape(example_watershed, res_ws1)
names(example_watershed_mod)

Note that a new column state appears in now in the sf object. We can check the effect by drawing the relative water content:

plot_variable(example_watershed_mod, variable = "soil_rwc_curr", r = r)

Now we can continue our simulation, in this case adding an extra month:

dates <- seq(as.Date("2001-02-01"), as.Date("2001-02-28"), by="day")
res_ws3 <- spwb_land(r, example_watershed_mod,
                    SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month",
                    watershed_control = ws_control, progress = FALSE)

The fact that no cell required initialization is an indication that we used an already initialized landscape.

Burn-in periods

Like other distributed hydrological models, watershed simulations with medfateland will normally require a burn-in period to allow soil moisture and aquifer levels to reach a dynamic equilibrium. We recommend users to use at least one or two years of burn-in period, but this will depend on the size of the watershed. In medfate we provide users with a copy of the example watershed, where burn-in period has already been simulated. This can be seen by inspecting the aquifer level:

data("example_watershed_burnin")
plot_variable(example_watershed_burnin, variable = "aquifer", r = r)

If we run a one-month simulation on this data set we can then compare the output before and after the burn-in period to illustrate its importance:

dates <- seq(as.Date("2001-01-01"), as.Date("2001-01-31"), by="day")
res_ws3 <- spwb_land(r, example_watershed_burnin,
                    SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month",
                    watershed_control = ws_control, progress = FALSE)
data.frame("before" = res_ws1$watershed_balance$WatershedExport, 
           "after" = res_ws3$watershed_balance$WatershedExport)

Simulations of watershed forest dynamics

Running growth_land() is very similar to running spwb_land(). However, a few things change when we want to simulate forest dynamics using fordyn_land(). Regarding the sf input, an additional column management_arguments may be defined to specify the forest management arguments (i.e. silviculture) of cells. Furthermore, the function does not allow choosing the temporal scale of summaries. Strong simplification of forest structure to dominant species will not normally make sense in this kind of simulation, since the focus is on forest dynamics.

A call to fordyn_land() for a single year is given here, as an example, starting from the initial example watershed:

res_ws4 <- fordyn_land(r, example_watershed,
                       SpParamsMED, examplemeteo,
                       watershed_control = ws_control, progress = FALSE)

Simulations using weather interpolation

Large watersheds will have spatial differences in climatic conditions like temperature, precipitation. Specifying a single weather data frame for all the watershed may be not suitable in this case. Specifying a different weather data frame for each watershed cell can also be a problem, if spatial resolution is high, due to the huge data requirements. A solution for this can be using interpolation on the fly, inside watershed simulations. This can be done by supplying an interpolator object (or a list of them), as defined in package meteoland. Here we use the example data provided in the package:

interpolator <- meteoland::with_meteo(meteoland_meteo_example, verbose = FALSE) |>
    meteoland::create_meteo_interpolator(params = defaultInterpolationParams())

Once we have this object, using it is straightforward:

res_ws5 <- spwb_land(r, example_watershed_burnin, SpParamsMED, 
                     meteo = interpolator, summary_frequency = "month",
                     watershed_control = ws_control, progress = FALSE)

Note that we did not define dates, which are taken from the interpolator data. If we plot the minimum temperature, we will appreciate the spatial variation in climate:

plot_summary(res_ws5$sf, variable = "MinTemperature", date = "2022-04-01", r = r)

For large watersheds and fine spatial resolution interpolation can become slow. One can then specify that interpolation is performed on a coarser grid, by using a watershed control parameter, for example:

ws_control$weather_aggregation_factor <- 3

To illustrate its effect, we repeat the previous simulation and plot the minimum temperature:

res_ws6 <- spwb_land(r, example_watershed_burnin, SpParamsMED, 
                     meteo = interpolator, summary_frequency = "month",
                     watershed_control = ws_control, progress = FALSE)
plot_summary(res_ws6$sf, variable = "MinTemperature", date = "2022-04-01", r = r)

Parallel simulations using subwatersheds

Simulations can be rather slow even for moderately-sized watersheds. For these reason, medfateland now incorporates the possibility to perform parallel simulations in subwatersheds, and then aggregate the results. To illustrate this feature we will use the data set of Bianya watershed (see Preparing inputs II: arbitrary locations; you can find the dataset in the medfateland GitHub repository).

We begin by loading the raster and sf inputs for Bianya:

r <- terra::rast("../intro/bianya_raster.tif")
sf <- readRDS("../intro/bianya.rds")

If we draw the elevation map, we will visually identify two subwatersheds:

plot_variable(sf, variable = "elevation", r = r)

Package medfateland includes function overland_routing() to statically illustrate how overland runoff processes are dealt with (i.e. distribution of runoff among neighbors, channel routing, etc.).

or <- overland_routing(r, sf)
head(or)

In this function we can specifically ask for sub-watersheds, as follows:

or <- overland_routing(r, sf, subwatersheds = TRUE,
                       max_overlap = 0.3)

In short, sub-watershed definition consists in: (a) finding the drainage basin of each outlet or channel cell; (b) aggregating drainage basins until the overlap is less than the specified parameter; and (c) deciding to which sub-watershed each border cell belongs to. This function identified three sub-watersheds, one of them being an isolated channel cell:

plot(or[, "subwatershed"])

Let's now illustrate how to perform watershed simulations with parallelization and subwatersheds. We begin by initializing the input (here we used a "buckets" soil hydrology to speed up calculations, but "single" would be more appropriate).

sf_init <- initialize_landscape(sf, SpParams = SpParamsMED,
                                local_control = defaultControl(soilDomains = "buckets"),
                                progress = FALSE)

Subwatershed definition is controlled via watershed control options as follows (simulation functions internally call overland_routing()):

ws_control <- default_watershed_control("tetis")
ws_control$tetis_parameters$subwatersheds <- TRUE
ws_control$tetis_parameters$max_overlap <- 0.3

Now, we are ready to launch the watershed simulation with parallelization. This consists in performing simulations for each subwatershed independently, aggregating the results and, finally, performing channel routing.

For simplicity, we only simulate five days. We ask for console output to see what the model is doing:

dates <- seq(as.Date("2022-04-01"), as.Date("2022-04-05"), by="day")
res_ws7 <- spwb_land(r, sf_init,
                    SpParamsMED, interpolator, dates = dates,
                    summary_frequency = "day", summary_blocks = "WaterBalance",
                    watershed_control = ws_control, progress = TRUE,
                    parallelize = TRUE)

As an example of the output, we show a map of woody plant transpiration (note that we asked for water balance components using summary_blocks = "WaterBalance":

plot_summary(res_ws7, variable = "Transpiration", date = "2022-04-01", r = r)

emf-creaf/medfateland documentation built on April 17, 2025, 5:43 a.m.