Watershed simulations
Miquel De Caceres
2025-04-08
Source:vignettes/runmodels/WatershedSimulations.Rmd
WatershedSimulations.RmdAim
The aim of this vignette is to illustrate how to use
medfateland (v. 2.6.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 (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## 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.75Channel 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## 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(soilDomains = "buckets"))## ℹ Creating 65 state objects for model 'spwb'.## ✔ Creating 65 state objects for model 'spwb'. [15ms]## ## • Transpiration mode [Granier: 65, Sperry: 0, Sureau: 0]## • Soil doimains [buckets: 65, single: 0, dual: 0]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]## ## • Transpiration mode [Granier: 65, Sperry: 0, Sureau: 0]## • Soil doimains [buckets: 65, single: 0, dual: 0]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 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)## [1] "watershed_control" "sf" "watershed_balance"
## [4] "watershed_soil_balance" "channel_export_m3s" "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 7 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 × 8
## geometry state aquifer snowpack summary result outlet
## * <POINT [m]> <list> <dbl> <dbl> <list> <list> <lgl>
## 1 (402630 4672570) <spwbInpt [19]> 123. 3.56 <dbl[…]> <NULL> FALSE
## 2 (402330 4672470) <aspwbInp [4]> 0.308 3.56 <dbl[…]> <NULL> FALSE
## 3 (402430 4672470) <spwbInpt [19]> 2.19 3.56 <dbl[…]> <NULL> FALSE
## 4 (402530 4672470) <spwbInpt [19]> 9.65 2.54 <dbl[…]> <NULL> FALSE
## 5 (402630 4672470) <spwbInpt [19]> 148. 2.57 <dbl[…]> <NULL> FALSE
## 6 (402730 4672470) <aspwbInp [4]> 868. 3.56 <dbl[…]> <NULL> TRUE
## 7 (402830 4672470) <aspwbInp [4]> 405. 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.03 3.56 <dbl[…]> <NULL> FALSE
## # ℹ 56 more rows
## # ℹ 1 more variable: outlet_backlog <list>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)## 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.00000 0.000000
## 2 1.889229 0.00000000 0.000000 0.00000 0.000000
## 3 0.000000 0.00000000 0.000000 0.00000 0.000000
## 4 5.024594 0.00000000 5.399648 53.29682 5.399648
## 5 1.403519 0.00000000 8.620277 63.95029 8.620277
## 6 12.447495 0.05090607 14.321824 68.00257 14.372730
## DeepDrainage CapillarityRise DeepAquiferLoss SoilEvaporation Transpiration
## 1 2.870064 0 0 0.3874890 0.2846843
## 2 1.766980 0 0 0.4675914 0.5045425
## 3 1.910902 0 0 0.3482580 0.4236354
## 4 2.275586 0 0 0.2006258 0.1903985
## 5 2.175750 0 0 0.2893012 0.5029694
## 6 2.335915 0 0 0.1855159 0.3603860
## HerbTranspiration InterflowBalance BaseflowBalance AquiferExfiltration
## 1 0 0.000000e+00 0.000000e+00 0
## 2 0 1.354136e-16 2.439126e-17 0
## 3 0 4.735269e-16 -5.298792e-17 0
## 4 0 1.703183e-16 2.775558e-17 0
## 5 0 1.350771e-15 -2.102695e-17 0
## 6 0 -2.974052e-15 7.611756e-17 0
## ChannelExport WatershedExport
## 1 0 0.000000
## 2 0 0.000000
## 3 0 0.000000
## 4 0 5.399648
## 5 0 8.620277
## 6 0 14.372730Values 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)## 6
## 2001-01-01 0.00000000
## 2001-01-02 0.00000000
## 2001-01-03 0.00000000
## 2001-01-04 0.04127060
## 2001-01-05 0.06588651
## 2001-01-06 0.10985367Accessing 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 SWE
## 2001-01-01 -3.203556 2.427977 31.14151 58.09884 16.65065 1.699603
## RWC SoilVol WTD DTA
## 2001-01-01 97.68028 525.6724 NA 14.33906Additional 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)] <- 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"
## [16] "outlet_backlog"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.09650658
## 3 0.00000000 0.09592942
## 4 5.39964799 0.09555595
## 5 8.62027675 0.09459706
## 6 14.37273028 0.12787286
## 7 9.72905753 0.09423911
## 8 17.11843924 0.09402956
## 9 12.80160239 0.09503276
## 10 5.84983613 0.10051036
## 11 6.84532385 0.10256741
## 12 10.27732596 0.10906450
## 13 7.93248397 0.11632446
## 14 7.90716679 0.12064020
## 15 5.75924214 0.12002444
## 16 3.96339743 3.04500877
## 17 3.39266678 3.00673669
## 18 2.74153725 2.33180505
## 19 1.97001006 0.82093735
## 20 4.72546190 0.22336146
## 21 1.59775872 0.26273819
## 22 0.79103009 0.22293555
## 23 0.33483610 0.17877538
## 24 0.30859359 0.17760500
## 25 0.28805641 0.17680538
## 26 0.24383072 0.17536188
## 27 0.22968842 0.17350457
## 28 0.20651256 0.17141103
## 29 0.18311761 0.16902006
## 30 0.08320080 0.16650391
## 31 0.03666704 0.16395753Simulations 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)
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:
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)## Simple feature collection with 6 features and 11 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 449799.3 ymin: 4678387 xmax: 450794.5 ymax: 4678387
## Projected CRS: ETRS89 / UTM zone 31N
## # A tibble: 6 × 12
## geometry elevation waterRank waterOrder queenNeigh waterQ
## <POINT [m]> <dbl> <int> <int> <list> <list>
## 1 (449799.3 4678387) 900. 401 2533 <int [4]> <dbl [4]>
## 2 (449998.3 4678387) 901. 399 2455 <int [5]> <dbl [5]>
## 3 (450197.4 4678387) 880. 446 2416 <int [5]> <dbl [5]>
## 4 (450396.4 4678387) 843. 554 2487 <int [5]> <dbl [5]>
## 5 (450595.5 4678387) 878. 454 2534 <int [5]> <dbl [5]>
## 6 (450794.5 4678387) 860. 508 2511 <int [4]> <dbl [4]>
## # ℹ 6 more variables: channel <lgl>, outlet <lgl>, target_outlet <int>,
## # distance_to_outlet <dbl>, outlet_backlog <list>, subwatersheds <lgl>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.3Now, 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)## ## ── Simulation of model 'spwb' over a watershed ─────────────────────────────────## ## ── INPUT CHECKING ──## ## ℹ Checking raster topology## ✔ Checking raster topology [13ms]## ## ℹ Checking 'sf' data columns## ℹ Column 'snowpack' was missing in 'sf'. Initializing empty snowpack.## ℹ Checking 'sf' data columnsℹ Minimum bedrock porosity set to 0.1%.
## ℹ Checking 'sf' data columnsℹ Column 'aquifer' was missing in 'sf'. Initializing empty aquifer.
## ℹ Checking 'sf' data columns✔ Checking 'sf' data columns [3s]
##
## ℹ Determining neighbors and overland routing for TETIS
## ✔ Determining neighbors and overland routing for TETIS [1.7s]
##
## • Hydrological model: TETIS
## • Number of grid cells: 3825 Number of target cells: 2573
## • Average cell area: 39575 m2, Total area: 15138 ha, Target area: 10183 ha
## • Cell land use [wildland: 2161 agriculture: 331 artificial: 78 rock: 0 water:
## 3]
## • Cells with soil: 2492
## • Number of days to simulate: 5
## • Number of temporal cell summaries: 5
## • Number of cells with daily model results requested: 0
## • Number of channel cells: 158
## • Number of outlet cells: 28
## • Number of subwatersheds: 3
## • Weather interpolation factor: 1
##
## ── INITIALISATION ──
##
## ℹ All state objects are already available for 'spwb'.
## • Transpiration mode [Granier: 2492, Sperry: 0, Sureau: 0]
## • Soil doimains [buckets: 2492, single: 0, dual: 0]
##
## ── PARALLEL SIMULATION of 3 SUB-WATERSHEDS in 3 NODES ──
##
## ── MERGING SUB-WATERSHED RESULTS ──
##
## ── CHANNEL ROUTING ──
##
## • Initial outlet backlog sum (m3): 0
## • Channel balance target (m3): 0 outlet change (m3): 0 backlog change (m3): 0
## • Final outlet backlog sum (m3): 0As 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)