Adding catchment data

The basis of the analysis is a rasterised Digital Elevation Model (DEM) of the catchment and a vectorised representation of the river network with attributes. Currently these can be in any format supported by the terra library.

However, within the calculations used for sink filling and flow routing the raster DEM is presumed to be projected so that is has square cells such that the difference between the cell centres (in meters) does not alter.

For Swindale the suitable DEM and channel files can be found using:

dem_file <- system.file("extdata/gis", "SwindaleDTM.tif", package="dynatop", mustWork = TRUE)
channel_file <- system.file("extdata/gis", "SwindaleRiverNetwork.gpkg", package="dynatop", mustWork = TRUE)

Alongside these a catchment outline must be provided to define the catchment boundary. In this example the catchment map is located at

catchment_file <- system.file("extdata/gis", "SwindaleBoundary.gpkg", package="dynatop", mustWork = TRUE)

To start the analysis the DEm and catchment boundary are initially added. The resolution and projection are taken from the DEM. In this case we add both as paths to the approiate files:

ctch$add_catchment(catchment_file,dem_file)

Adding river channel data is more complex. The add_channel method requires a SpatVector object (or a file name that can be loaded as a SpatVect object) containing lines defining the river centre. Each vector object is treated as a length of river channel which requires the following properties

• uid - a label for the channel length
• endNode - a label for the downstream end of the river length
• startNode - a label for the upstream end of the river length
• width - width of the channel
• depth - depth of the channel
• slope - bed slope of the channel

Since it is possible that these properties are present in a data file under different names some basic preprocessing may be required. The convert_channel function is designed to help with this. To illustrate this let us examine the river network for Swindale

chn <- terra::vect(channel_file)
head(chn)
#> class       : SpatVector
#> geometry    : lines
#> dimensions  : 6, 11  (geometries, attributes)
#> extent      : 349884.5, 351675.7, 508614.5, 513074.7  (xmin, xmax, ymin, ymax)
#> coord. ref. : OSGB36 / British National Grid (EPSG:27700)
#> names       :           name1      identifier       startNode         endNode        form         flow fictitious length name2 sinkdepth Shape_Leng
#> type        :           <chr>           <chr>           <chr>           <chr>       <chr>        <chr>      <chr>  <num> <chr>     <num>      <num>
#> values      :   Hawthorn Gill 16D0AC09-E0B6-~ 730F01D2-0F4F-~ 72B1A40B-E106-~ inlandRiver in direction      false    431    NA        -1    431.304
#>               Little Mosedal~ 643017BB-01BC-~ 9B06188C-F4C7-~ A19E0D3A-2FA1-~ inlandRiver in direction      false    513    NA        -1    513.128
#>                 Mosedale Beck D2C2703E-A32F-~ B03DB3BD-0E33-~ D649B60C-E631-~ inlandRiver in direction      false    740    NA        -1    739.764
#>               ...

Some of the main properties are present under appropriate names (startNode, endNode) but the remainder are missing. The suitability of a channel can be tested by

check_channel(chn)
#> Error in `check_channel()`:
#> ! uid property is missing

The convert_channelfunction addresses these shortcomings by renaming and checking the approriate type (e.g. character , double, integer) of the variabes. It does not currently set default values for the channle properties, except for slope. To illustrate add vector of missing values to the chn.

chn$slope <- NA
chn$slope[1] <- 0
chn$class <- 1
chn$width <- 1
chn$depth <- 0.1

then running convert_channel will run and repopulate the slope values with a default (1e-3)

chn <- convert_channel(chn,
                       property_names=c(uid = "identifier",
                                        startNode = "startNode",
                                        endNode = "endNode",
                                        slope = "slope",
                                        class = "class",
                                        width="width",
                                        depth="depth"))
#> Warning in convert_channel(chn, property_names = c(uid = "identifier",
#> startNode = "startNode", : None finite slope values replaced by minimum
#> Warning in convert_channel(chn, property_names = c(uid = "identifier",
#> startNode = "startNode", : Shallow slope values replaced by minimum

As shown above the convert_channel function takes a named vector giving the variable names to be use for the required channel properties. More advanced processof channel is available allowing them to be simplied by merging lines (simplify_channel) or be trimmed to certain outlets (trim_channel).

The river network can then be added to the project by

ctch$add_channel(chn)

Getting and plotting catchment information

The dynatopGIS class has methods for returning and plotting the GIS data in the project. The names of all the different GIS layers stored is returned by

ctch$get_layer()
#> [1] "catchment"     "dem"           "channel_fd"    "channel_slp"  
#> [5] "channel_class" "hru"           "channel_width" "channel_depth"

These can be plotted (with or without the channel), for example

ctch$plot_layer("dem", add_channel=TRUE)

or returned, for example

ctch$get_layer("dem")
#> class       : SpatRaster
#> size        : 163, 124, 1  (nrow, ncol, nlyr)
#> resolution  : 40, 40  (x, y)
#> extent      : 347734, 352694, 507244, 513764  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
#> source(s)   : memory
#> varname     : SwindaleDTM
#> name        :        dem
#> min value   : 262.800446
#> max value   : 710.753296

All layers are returned as SpatRast objects.

Filling sinks

For the hill slope to be connected to the river network all DEM cells must drain to those that intersect with the river network. The algorithm implemented in the sink_fill method ensures this is the case.

The algorithm of used by the sink_fill method is iterative and the execution time of the function is limited by capping the maximum number of iterations. If this limit is reached without completion the method can call again with the “hot start” option to continue from where it finished.

For Swindale, where the example DEM is already partially filled the algorithm only alters a small area near the foot of the catchment.

ctch$sink_fill()

terra::plot( ctch$get_layer('filled_dem') - ctch$get_layer('dem'),
            main="Changes to height")

Determining Ordering

The computational scheme in the dynatop package works with an ordered sequence of HRUs constructed such that the sequence moves downslope to catchment outlet. This is achieved by banding the channel reaches and hillslope cells such that the catchment outlet(s) are in band 0, those cells or reaches draining only into band 1 are in band 2 and so forth. banding is achieved by the following call

ctch$order()
ctch$plot_layer("band")

Adding additional layers

Properties may come in additional GIS layers. To demonstrate the addition of an additional layer we will extract the filled dem

tmp <- ctch$get_layer("filled_dem")

then separate it into 6 layers representing different altitude bands.

s <- terra::classify(tmp,5)
names(s) <- "alt_band"

The resulting raster object can now be added to the project with

ctch$add_layer(s)
ctch$get_layer()
#>  [1] "catchment"     "dem"           "channel_fd"    "channel_slp"  
#>  [5] "channel_class" "hru"           "channel_width" "channel_depth"
#>  [9] "filled_dem"    "band"          "alt_band"

Note the use of the name to identify the layer.

Producing a catchment classification

A dynatop model uses a catchment classification to assign properties to each hydrological response unit (HRU). In this example we will use a single layer, the alt_band layer just computed, to provide the classification. The classification can be produced by

ctch$combine_layers("class_alt_band","alt_band")

In producing the classification a json fiel containing the properties from each combined layers associated with the class is generated and saved in the project directory.

Making the model

Once a clasification is available the model is generatd by a single function call

mdlRoot <- file.path( demo_dir, "test_mdl")
ctch$create_model(mdlRoot,class_layer="class_alt_band")

Two files are produced.

list.files(demo_dir, "test_mdl")
#> [1] "test_mdl_param.json" "test_mdl.gpkg"

The first file test_mdl.gpkg file contains the spatial location and physical properties and class of each HRU. The second file test_mdl.json contains example parameters for each class.

Preparing input data

The input to the model is expected to take the form of an data.frame with a column called dateTime give the times of the observed data using a constant time step. The remaining columns contain data series whose column names are found in the ‘precip’ and ‘pet’ variables of the HRUs descriptions in the model.

The discharge, precipitation and potential evapotranspiration (PET) inputs for Swindale are available on a 15 minute time step.

obs_file <- system.file("extdata","mdl","SwindaleObs.csv",package="dynatop",mustWork=TRUE)
obs <- read_obs(obs_file)
head(obs)
#>              dateTime flow precip          pet
#> 1 2009-11-18 16:00:00 2.78  4e-04 8.878467e-06
#> 2 2009-11-18 16:15:00 2.80  2e-04 4.762229e-06
#> 3 2009-11-18 16:30:00 2.85  4e-04 8.585708e-07
#> 4 2009-11-18 16:45:00 2.94  2e-04 0.000000e+00
#> 5 2009-11-18 17:00:00 3.02  2e-04 0.000000e+00
#> 6 2009-11-18 17:15:00 3.15  2e-04 0.000000e+00

Note the discharge is in m\(^{3}\)/s while the precipitation and PET are in m accumulated since the previous observation time.

To use the data with the model we need to set the names of the time series inputs within the model. In this case this is already done as can be seen by inspecting the precip and pet values for the individual HRU’s.

geo <- read_geo( file.path(demo_dir,"test_mdl.gpkg") ) ## read in model
unique(geo$properties$precip) ## precipitation time series names
#> NULL
unique(geo$properties$pet) ## PET time series names
#> NULL

Altering parameters

The parameter values are stored within the parameter file. Altering parameter values requires changing their values, either in the file, or in the object once read in. For this example all HRUs, with the exception of the channel HRUs will have the same parameter values. Changing the parameters can be carried out using a function.

param <- read_param( file.path(demo_dir,"test_mdl_param.json") )
fchange <- function(x){
    x$K_s <- 9/3600
    m <- 1
    if(!is.na(x$channel_class) & x$channel_class == 1){
        ## river
        x$K_s <- 0 ##  force water to surface
    }
    x$r_uz_sz_0 <- 1.755582e-07 ## initialisation parameter - the initial outflow divided by catchment area
    return(x)
}
p <- lapply(param,fchange)

Creating the dynatop Object

Simulations are performed by embedding the model and the observed data into a dynatop object. First the object is created using the model and parameters

ctch_mdl <- dynatop$new(geo,p)

This step performs some basic checks on the model for consistency.

Defining outputs

To run the model requires not both the observed data and a definition of the output variables to collect.

The output definition is adata.frame with four columns

defn <- data.frame(name = "gauge_flow",
                   hru = 1,
                   scale = 1,
                   flux = "q_sf")

Initialising the states

The initial states of the simulations can be specified directly or estimated using a quazi steady state approach based on an initial recharge rate (the parameter r_us_sz_0). To use the latter

ctch_mdl$initialise()

Simulation

The simulation can now be performed by passing the observed data and output definition to the model.

sim1 <- ctch_mdl$sim(obs,defn)
#> Number of threads 1

Note that the states of the system are now those at the end of the simulation for example:

head(ctch_mdl$get_states())
#>   hru        s_sf s_rz s_uz s_sz
#> 1   0 3464.269455  160    0    0
#> 2   1 3453.029462  160    0    0
#> 3   2    8.500562  160    0    0
#> 4   3 3441.263910  160    0    0
#> 5   4 3429.424238  160    0    0
#> 6   5    3.590769  160    0    0

Rerunning the simulation with the new initial states will of course produce different results. Output for the above examples can be plotted against observed discharge for comparison as follows:

sim2 <- ctch_mdl$sim(obs,defn)
#> Number of threads 1
out <- Reduce( cbind, list("obs" = obs$flow, "sim1" = sim1$gauge_flow, "sim2"=sim2$gauge_flow))
matplot(obs$dateTime,out)