Watershed simulations
Miquel De Caceres
2024-11-08
Source:vignettes/runmodels/WatershedSimulations.Rmd
WatershedSimulations.RmdAim
The aim of this vignette is to illustrate how to use
medfateland (v. 2.5.0) 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.
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## Simple feature collection with 66 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 401430 ymin: 4671870 xmax: 402830 ymax: 4672570
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 15
## geometry id elevation slope aspect land_cover_type
## * <POINT [m]> <int> <dbl> <dbl> <dbl> <chr>
## 1 (402630 4672570) 1 1162 11.3 79.2 wildland
## 2 (402330 4672470) 2 1214 12.4 98.7 agriculture
## 3 (402430 4672470) 3 1197 10.4 102. wildland
## 4 (402530 4672470) 4 1180 8.12 83.3 wildland
## 5 (402630 4672470) 5 1164 13.9 96.8 wildland
## 6 (402730 4672470) 6 1146 11.2 8.47 agriculture
## 7 (402830 4672470) 7 1153 9.26 356. agriculture
## 8 (402230 4672370) 8 1237 14.5 75.1 wildland
## 9 (402330 4672370) 9 1213 13.2 78.7 wildland
## 10 (402430 4672370) 10 1198 8.56 75.6 agriculture
## # ℹ 56 more rows
## # ℹ 9 more variables: forest <list>, soil <list>, state <list>,
## # depth_to_bedrock <dbl>, bedrock_conductivity <dbl>, bedrock_porosity <dbl>,
## # snowpack <dbl>, aquifer <dbl>, crop_factor <dbl>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)##
## agriculture rock wildland
## 17 1 48Local 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.75Grid 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## class : SpatRaster
## dimensions : 8, 15, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 401380, 402880, 4671820, 4672620 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 31N (EPSG:32631)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())## ℹ Creating 65 state objects for model 'spwb'.## ✔ Creating 65 state objects for model 'spwb'. [18ms]##
example_init## Simple feature collection with 66 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 401430 ymin: 4671870 xmax: 402830 ymax: 4672570
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 15
## geometry id elevation slope aspect land_cover_type
## * <POINT [m]> <int> <dbl> <dbl> <dbl> <chr>
## 1 (402630 4672570) 1 1162 11.3 79.2 wildland
## 2 (402330 4672470) 2 1214 12.4 98.7 agriculture
## 3 (402430 4672470) 3 1197 10.4 102. wildland
## 4 (402530 4672470) 4 1180 8.12 83.3 wildland
## 5 (402630 4672470) 5 1164 13.9 96.8 wildland
## 6 (402730 4672470) 6 1146 11.2 8.47 agriculture
## 7 (402830 4672470) 7 1153 9.26 356. agriculture
## 8 (402230 4672370) 8 1237 14.5 75.1 wildland
## 9 (402330 4672370) 9 1213 13.2 78.7 wildland
## 10 (402430 4672370) 10 1198 8.56 75.6 agriculture
## # ℹ 56 more rows
## # ℹ 9 more variables: forest <list>, soil <list>, state <list>,
## # depth_to_bedrock <dbl>, bedrock_conductivity <dbl>, bedrock_porosity <dbl>,
## # snowpack <dbl>, aquifer <dbl>, crop_factor <dbl>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)## ℹ Creating 65 state objects for model 'spwb'.## ✔ Creating 65 state objects for model 'spwb'. [5ms]##
example_simplified## Simple feature collection with 66 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 401430 ymin: 4671870 xmax: 402830 ymax: 4672570
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 15
## geometry id elevation slope aspect land_cover_type
## * <POINT [m]> <int> <dbl> <dbl> <dbl> <chr>
## 1 (402630 4672570) 1 1162 11.3 79.2 wildland
## 2 (402330 4672470) 2 1214 12.4 98.7 agriculture
## 3 (402430 4672470) 3 1197 10.4 102. wildland
## 4 (402530 4672470) 4 1180 8.12 83.3 wildland
## 5 (402630 4672470) 5 1164 13.9 96.8 wildland
## 6 (402730 4672470) 6 1146 11.2 8.47 agriculture
## 7 (402830 4672470) 7 1153 9.26 356. agriculture
## 8 (402230 4672370) 8 1237 14.5 75.1 wildland
## 9 (402330 4672370) 9 1213 13.2 78.7 wildland
## 10 (402430 4672370) 10 1198 8.56 75.6 agriculture
## # ℹ 56 more rows
## # ℹ 9 more variables: forest <list>, soil <list>, state <list>,
## # depth_to_bedrock <dbl>, bedrock_conductivity <dbl>, bedrock_porosity <dbl>,
## # snowpack <dbl>, aquifer <dbl>, crop_factor <dbl>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 results at
coarser temporal scales, to reduce the amount of memory in spatial
results.
Structure of simulation outputs
Function spwb_land() and growth_land()
return a list with the following elements:
names(res_ws1)## [1] "watershed_control" "sf" "watershed_balance"
## [4] "watershed_soil_balance" "outlet_export_m3s"Where sf is an object of class sf,
analogous to those of functions *_spatial():
res_ws1$sf## Simple feature collection with 66 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 401430 ymin: 4671870 xmax: 402830 ymax: 4672570
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 7
## geometry state aquifer snowpack summary result outlet
## <POINT [m]> <list> <dbl> <dbl> <list> <list> <lgl>
## 1 (402630 4672570) <spwbInpt [19]> 127. 3.56 <dbl[…]> <NULL> FALSE
## 2 (402330 4672470) <aspwbInp [4]> 0.308 3.56 <dbl[…]> <NULL> FALSE
## 3 (402430 4672470) <spwbInpt [19]> 2.21 3.56 <dbl[…]> <NULL> FALSE
## 4 (402530 4672470) <spwbInpt [19]> 10.1 2.54 <dbl[…]> <NULL> FALSE
## 5 (402630 4672470) <spwbInpt [19]> 155. 2.57 <dbl[…]> <NULL> FALSE
## 6 (402730 4672470) <aspwbInp [4]> 876. 3.56 <dbl[…]> <NULL> TRUE
## 7 (402830 4672470) <aspwbInp [4]> 410. 3.56 <dbl[…]> <NULL> FALSE
## 8 (402230 4672370) <spwbInpt [19]> 0.261 2.53 <dbl[…]> <NULL> FALSE
## 9 (402330 4672370) <spwbInpt [19]> 1.45 2.81 <dbl[…]> <NULL> FALSE
## 10 (402430 4672370) <aspwbInp [4]> 6.11 3.56 <dbl[…]> <NULL> FALSE
## # ℹ 56 more rowsColumns state, aquifer and
snowpack contain state variables, whereas
summary contains temporal water balance 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)## dates Precipitation Rain Snow Snowmelt Interception NetRain
## 1 2001-01-01 4.869109 4.869109 0 0 0.7754218 4.093687
## 2 2001-01-02 2.498292 2.498292 0 0 0.6090627 1.889229
## 3 2001-01-03 0.000000 0.000000 0 0 0.0000000 0.000000
## 4 2001-01-04 5.796973 5.796973 0 0 0.7723790 5.024594
## 5 2001-01-05 1.884401 1.884401 0 0 0.4808826 1.403519
## 6 2001-01-06 13.359801 13.359801 0 0 0.8613997 12.498401
## Infiltration InfiltrationExcess SaturationExcess CellRunon CellRunoff
## 1 4.093687 0.00000000 0.000000 0.00000000 0.000000
## 2 1.889229 0.00000000 0.000000 0.00000000 0.000000
## 3 0.000000 0.00000000 0.000000 0.00000000 0.000000
## 4 5.024594 0.00000000 5.399544 0.00000000 5.399544
## 5 1.403519 0.00000000 8.620346 0.00000000 8.620346
## 6 12.498401 0.05090607 14.332811 0.05090607 14.383717
## DeepDrainage CapillarityRise DeepAquiferLoss SoilEvaporation Transpiration
## 1 2.870060 0 0 0.3874890 0.2846970
## 2 1.766971 0 0 0.4675914 0.5045596
## 3 1.910884 0 0 0.3482575 0.4236524
## 4 2.275585 0 0 0.2006254 0.1904111
## 5 2.175753 0 0 0.2893019 0.5026577
## 6 2.335919 0 0 0.1855165 0.3597430
## HerbTranspiration InterflowBalance BaseflowBalance AquiferExfiltration
## 1 0 0.000000e+00 0.000000e+00 0
## 2 0 -1.859834e-16 3.364312e-17 0
## 3 0 -7.372049e-16 2.523234e-18 0
## 4 0 3.462298e-15 -3.953067e-17 0
## 5 0 -1.364229e-15 1.261617e-17 0
## 6 0 3.606543e-15 5.803439e-17 0
## WatershedExport
## 1 0.000000
## 2 0.000000
## 3 0.000000
## 4 5.399544
## 5 8.620346
## 6 14.383717Values 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,
outlet_export_m3 contains the volume reaching each outlet
cell per day:
head(res_ws1$outlet_export_m3)## 6
## 2001-01-01 0.00000000
## 2001-01-02 0.00000000
## 2001-01-03 0.00000000
## 2001-01-04 0.04126980
## 2001-01-05 0.06588704
## 2001-01-06 0.10993764Accessing 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]]## MinTemperature MaxTemperature PET Rain Snow Snowmelt
## 2001-01-01 -3.203556 2.427977 31.14151 58.09884 16.65065 13.09301
## Interception NetRain Infiltration InfiltrationExcess
## 2001-01-01 23.59187 34.50697 47.59998 0
## SaturationExcess Runon Runoff DeepDrainage CapillarityRise
## 2001-01-01 0 0 0 55.82266 0
## DeepAquiferLoss SoilEvaporation Transpiration HerbTranspiration
## 2001-01-01 0 2.792085 9.971227 0
## InterflowInput InterflowOutput InterflowBalance BaseflowInput
## 2001-01-01 1040.115 1077.014 -36.89932 246.1273
## BaseflowOutput BaseflowBalance AquiferExfiltration SWE SoilVol
## 2001-01-01 174.5932 71.53411 0 1.699603 530.3081
## RWC WTD DTA
## 2001-01-01 98.54167 NA 14.32596Several plots 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)] <- TRUEIf 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)## [1] "spwb" "list"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)## [1] "geometry" "id" "elevation"
## [4] "slope" "aspect" "land_cover_type"
## [7] "forest" "soil" "state"
## [10] "depth_to_bedrock" "bedrock_conductivity" "bedrock_porosity"
## [13] "snowpack" "aquifer" "crop_factor"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)## before after
## 1 0.00000000 0.09798985
## 2 0.00000000 0.09650657
## 3 0.00000000 0.09592939
## 4 5.39954370 0.09555589
## 5 8.62034607 0.09459699
## 6 14.38371718 0.12787276
## 7 9.72970344 0.09423899
## 8 17.12959233 0.09402957
## 9 12.80189935 0.09503280
## 10 5.86754069 0.10051038
## 11 6.85004331 0.10256739
## 12 10.04891979 0.10906381
## 13 7.92767125 0.11632381
## 14 7.90721531 0.12063961
## 15 5.75790391 0.12002380
## 16 3.90833523 3.04502260
## 17 3.01041130 3.00683826
## 18 2.89642554 2.33187488
## 19 2.09476284 0.82091581
## 20 4.78444693 0.22336267
## 21 1.66854021 0.26273885
## 22 0.76462819 0.22293633
## 23 0.35263182 0.17877625
## 24 0.44476490 0.17760578
## 25 0.28943777 0.17680596
## 26 0.24380897 0.17536251
## 27 0.22966862 0.17350546
## 28 0.20649455 0.17141156
## 29 0.18310116 0.16901969
## 30 0.08320216 0.16650239
## 31 0.03666808 0.16395547Simulations 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())## ℹ Creating interpolator...## • Calculating smoothed variables...## • Updating intial_Rp parameter with the actual stations mean distance...## ✔ Interpolator created.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 <- 3To 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)