Aim

The aim of this vignette is to illustrate how to use medfateland (v. 2.2.1) 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 13 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 × 14
##            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
## # ℹ 8 more variables: forest <list>, soil <list>, state <list>,
## #   depth_to_bedrock <dbl>, bedrock_conductivity <dbl>, bedrock_porosity <dbl>,
## #   snowpack <dbl>, aquifer <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          48

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

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")

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_watershed,
                    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] "sf"                     "watershed_balance"      "watershed_soil_balance"
## [4] "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 [17]> 0.0404      3.56 <dbl[…]> <NULL> FALSE 
##  2 (402330 4672470) <aspwbInp [3]>  0.198       3.56 <dbl[…]> <NULL> FALSE 
##  3 (402430 4672470) <spwbInpt [17]> 0.0566      3.56 <dbl[…]> <NULL> FALSE 
##  4 (402530 4672470) <spwbInpt [17]> 0.0482      2.56 <dbl[…]> <NULL> FALSE 
##  5 (402630 4672470) <spwbInpt [17]> 0.00994     2.57 <dbl[…]> <NULL> FALSE 
##  6 (402730 4672470) <aspwbInp [3]>  0.243       3.56 <dbl[…]> <NULL> TRUE  
##  7 (402830 4672470) <aspwbInp [3]>  0.165       3.56 <dbl[…]> <NULL> FALSE 
##  8 (402230 4672370) <spwbInpt [17]> 0.00275     2.84 <dbl[…]> <NULL> FALSE 
##  9 (402330 4672370) <spwbInpt [17]> 0.00419     2.97 <dbl[…]> <NULL> FALSE 
## 10 (402430 4672370) <aspwbInp [3]>  0.239       3.56 <dbl[…]> <NULL> FALSE 
## # ℹ 56 more rows

Columns 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.7900101  4.079099
## 2 2001-01-02      2.498292  2.498292    0        0    0.6919287  1.806363
## 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.7855456  5.011427
## 5 2001-01-05      1.884401  1.884401    0        0    0.5571451  1.327256
## 6 2001-01-06     13.359801 13.359801    0        0    0.8937189 12.466082
##   Infiltration InfiltrationExcess SaturationExcess     Runoff DeepDrainage
## 1     4.005324                  0                0 0.07377438  0.003545635
## 2     1.768510                  0                0 0.03785290  0.003504170
## 3     0.000000                  0                0 0.00000000  0.003465461
## 4     4.923594                  0                0 0.08783292  0.003787363
## 5     1.298705                  0                0 0.02855154  0.003944415
## 6    12.263661                  0                0 0.20242122  0.004965145
##   CapillarityRise SoilEvaporation Transpiration HerbTranspiration
## 1    2.247422e-05       0.3867130     0.2793608                 0
## 2    1.757712e-04       0.1356762     0.4955779                 0
## 3    5.613658e-04       0.1844224     0.4164836                 0
## 4    2.330954e-05       0.1257523     0.1874572                 0
## 5    7.269397e-04       0.1448613     0.5090040                 0
## 6    4.064161e-04       0.1293526     0.3742005                 0
##   InterflowBalance BaseflowBalance AquiferExfiltration WatershedExport
## 1     0.000000e+00    0.000000e+00                   0      0.07377438
## 2    -1.443806e-23   -1.673341e-31                   0      0.03785290
## 3    -1.604229e-23   -2.581727e-30                   0      0.00000000
## 4    -1.411722e-22    9.561950e-30                   0      0.08783292
## 5    -5.903563e-22    5.507683e-29                   0      0.02855154
## 6    -6.416916e-22   -1.223930e-29                   0      0.20242122

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, outlet_export_m3 contains the volume reaching each outlet cell per day:

head(res_ws1$outlet_export_m3)
##                       6
## 2001-01-01 0.0005638724
## 2001-01-02 0.0002893174
## 2001-01-03 0.0000000000
## 2001-01-04 0.0006713247
## 2001-01-05 0.0002182251
## 2001-01-06 0.0015471461

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]]
##            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      17.1751 40.92374     54.01675                  0
##            SaturationExcess Runoff DeepDrainage CapillarityRise SoilEvaporation
## 2001-01-01                0      0     0.066837      0.02646211        1.693804
##            Transpiration HerbTranspiration InterflowInput InterflowOutput
## 2001-01-01      6.765493                 0     0.06753704      0.03424541
##            InterflowBalance BaseflowInput BaseflowOutput BaseflowBalance
## 2001-01-01       0.03329163  2.804731e-10    1.50188e-10    1.302852e-10
##            AquiferExfiltration      SWE  SoilVol      RWC      WTD     DTA
## 2001-01-01                   0 1.649861 569.2058 105.7696 3625.407 15.4995

Several 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_watershed$result_cell <- FALSE
example_watershed$result_cell[c(3,9)] <- TRUE

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

res_ws1 <- spwb_land(r, example_watershed,
                    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] "result_cell"

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 = "soilrwc", 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.073774376 0.3731846
## 2  0.037852903 0.3327946
## 3  0.000000000 0.2715065
## 4  0.087832923 0.3926512
## 5  0.028551537 0.3233051
## 6  0.202421222 0.5931811
## 7  0.000000000 0.3305984
## 8  0.000000000 0.2612705
## 9  0.000000000 0.2580586
## 10 0.141107695 0.4860296
## 11 0.090092361 0.4547424
## 12 0.102705122 0.4643604
## 13 0.001455222 0.3007647
## 14 0.000000000 0.2752761
## 15 0.000000000 0.2753467
## 16 0.000000000 0.2716527
## 17 0.000000000 0.2674355
## 18 0.000000000 0.2645210
## 19 0.000000000 0.2637802
## 20 0.000000000 0.2644819
## 21 0.122633540 0.4753295
## 22 0.136742721 0.5118090
## 23 0.007795713 0.3039228
## 24 0.017187861 0.2902332
## 25 0.018558367 0.2997559
## 26 0.000000000 0.2814159
## 27 0.009952749 0.2797189
## 28 0.000000000 0.2753959
## 29 0.000000000 0.2738947
## 30 0.000000000 0.2687231
## 31 0.000000000 0.2587818

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. A call to fordyn_land() for a single year is given here starting from a watershed after burnin:

res_ws4 <- fordyn_land(r, example_watershed_burnin,
                       SpParamsMED, examplemeteo,
                       watershed_control = ws_control, progress = TRUE)
## 
## ── Simulation of model 'fordyn' over a watershed ───────────────────────────────
##  Checking topology
##  Checking topology [11ms]
## 
##  Checking 'sf' data
## • Hydrological model: TETIS
##  Checking 'sf' data• Number of grid cells: 120 Number of target cells: 66
##  Checking 'sf' data• Average cell area: 10006 m2, Total area: 120 ha, Target area: 66 ha
##  Checking 'sf' data• Cell land use wildland: 48 agriculture: 17 artificial: 0 rock: 1 water: 0
##  Checking 'sf' data• Cells with soil: 65
##  Checking 'sf' data• Number of years to simulate: 1
##  Checking 'sf' data
##  Checking 'sf' data── Initialisation 
##  Checking 'sf' data Creating 65 input objects for model 'growth'.
##  Creating 65 input objects for model 'growth'. [5ms]
## 
##  Checking 'sf' dataStands ■■■■■■■■■■■■■■■■■■■■■■■           74% | ETA:  1s
## Stands ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■   98% | ETA:  0s
##  Checking 'sf' data Checking 'sf' data [3.4s]
## 
##  Initializing 'fordyn' output tables
##  Initializing 'fordyn' output tables [103ms]
## 
## 
## ── Simulating year 2001 [1/1] 
## • Growth/mortality
## Daily simulations    1% | ETA:  4m
## Daily simulations ■■                                 2% | ETA:  4m
## Daily simulations ■■                                 3% | ETA:  4m
## Daily simulations ■■                                 5% | ETA:  4m
## Daily simulations ■■■                                6% | ETA:  3m
## Daily simulations ■■■                                8% | ETA:  3m
## Daily simulations ■■■■                               9% | ETA:  3m
## Daily simulations ■■■■                              11% | ETA:  3m
## Daily simulations ■■■■■                             12% | ETA:  3m
## Daily simulations ■■■■■                             13% | ETA:  3m
## Daily simulations ■■■■■                             15% | ETA:  3m
## Daily simulations ■■■■■■                            16% | ETA:  3m
## Daily simulations ■■■■■■                            18% | ETA:  3m
## Daily simulations ■■■■■■■                           19% | ETA:  3m
## Daily simulations ■■■■■■■                           20% | ETA:  3m
## Daily simulations ■■■■■■■                           21% | ETA:  3m
## Daily simulations ■■■■■■■■                          22% | ETA:  3m
## Daily simulations ■■■■■■■■                          24% | ETA:  3m
## Daily simulations ■■■■■■■■                          25% | ETA:  3m
## Daily simulations ■■■■■■■■■                         26% | ETA:  3m
## Daily simulations ■■■■■■■■■                         27% | ETA:  3m
## Daily simulations ■■■■■■■■■                         28% | ETA:  3m
## Daily simulations ■■■■■■■■■■                        29% | ETA:  3m
## Daily simulations ■■■■■■■■■■                        30% | ETA:  3m
## Daily simulations ■■■■■■■■■■                        32% | ETA:  3m
## Daily simulations ■■■■■■■■■■■                       33% | ETA:  3m
## Daily simulations ■■■■■■■■■■■                       34% | ETA:  3m
## Daily simulations ■■■■■■■■■■■■                      35% | ETA:  3m
## Daily simulations ■■■■■■■■■■■■                      36% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■                      38% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■                     39% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■                     41% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■                    42% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■                    43% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■                    45% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■                   46% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■                   47% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■                  49% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■                  50% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■                 52% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■                 53% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■                 54% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■■                56% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■■                57% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■■■               58% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■■■               60% | ETA:  2m
## Daily simulations ■■■■■■■■■■■■■■■■■■■               61% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■              63% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■              64% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■             66% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■             67% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■            70% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■           72% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■           73% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■           75% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■          76% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■          77% | ETA:  1m
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■         79% | ETA: 47s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■         80% | ETA: 44s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■         82% | ETA: 41s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■        83% | ETA: 37s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■        84% | ETA: 35s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■       86% | ETA: 32s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■       87% | ETA: 29s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■       88% | ETA: 27s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      89% | ETA: 25s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      90% | ETA: 22s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■      92% | ETA: 19s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     93% | ETA: 16s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■     94% | ETA: 14s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    95% | ETA: 10s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    97% | ETA:  7s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■    98% | ETA:  4s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■   99% | ETA:  1s
## Daily simulations ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s
## • Seed production/dispersal
## • Management/recruitment/resprouting

In this case, the sf part of the output contains additional columns, analogous to those of fordyn_scenario().

names(res_ws4$sf)
##  [1] "geometry"         "state"            "aquifer"          "snowpack"        
##  [5] "summary"          "forest"           "tree_table"       "shrub_table"     
##  [9] "dead_tree_table"  "dead_shrub_table" "cut_tree_table"   "cut_shrub_table"

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 <- 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)