1 Introduction

This worked example applies the dynatop modeling software to the River Eden catchment in Cumbria. It presumes that you have at least a basic knowledge of hydrology and ideally some experience of R.

Useful primers for R can be found on multiple websites such as that of Hadley Wickham or on CRAN but complete knowledge of their contents is not expected.

1.1 Prerequisits

To follow the worked example the course is based requires

1.1.1 An installation of R 4.6 or above

Installation instructions for R are available from CRAN for Linux, Windows and Mac. Full details of the R version and packages used in building this version of the training material can be found here.

1.1.2 A suitable text editor for altering the R scripts.

While it is possible to complete the course using the built in R GUI’s for Windows and Mac it is recommended to install an editor with syntax highlighting. Popular choices include the Rstudio IDE, VSCode (Instructions here) or Emacs with ESS all of which provide an integrated development environment.

1.1.3 An installed copy of the dynatop package

The dynatop package and any dependencies can be install in the standard way for R packages. At the R command prompt type

install.packages('dynatop', repos =c('https://waternumbers.r-universe.dev','https://cloud.r-project.org'))

1.1.4 Data and Scripts

Download the data for the examples as a zip file. Extracting this should give a directory “eden_data” with four subdirectories. In the code provided is is presumed that the working directory of the R session is “eden_data”.

To save copying and pasting from the web pages the R code used in the training course is available. These scripts should be used alongside the course material.

2 River Eden

The training course works through an example application based on the Eden Catchment in Cumbria, UK (see Figure below). At the Sheepmount gauging station downstream of Carlisle, which is the lowest monitored point in the catchment, the total area is 2400 km\(^2\).

The River Eden and its main tributaries including the Rivers Eamont, Irthing, Petteril and Caldew link a number of significant areas for runoff generation which include the northern and eastern parts of the Lake District; areas of the Northern Pennines and Kielder Forest. These areas are typified by high annual rainfall totals and steep terrain. Moreover the confluences of two of these main tributaries lie within the Carlisle urban area resulting in a complex situation where significant inflows on any tributary may result in flooding. In January 2005 Carlisle was flooded by an event with an estimated 0.0067 Annual Exceedance Probability (AEP) at the Sheepmount Gauge (Clarke 2005; Day 2005). Subsequently to this events in 2010, and particularly 2009, came close to over topping the main flood defences. The significance of the 2005 flood, combined with a post event survey of 263 water or wrack marks, has resulted in multiple inundation model studies (Neal et al. 2009; Horritt et al. 2011; Fewtrell et al. 2011) and modelling (Leedal et al. 2013, Smith et al. 2013).

Map of the River Eden Catchment

Figure 2.1: Map of the River Eden Catchment

2.1 Available Data

The directory “eden_data/unprocessed” contains the starting data given in the table below

File Contents
Eden_Flow.csv Discharge data (m\(^3\)/s) for the Sheepmount Gauge provided by the environment Agency
Eden_Gauge_Sites.csv Locations of selected gauges in the Eden Catchment from the NRFA
Eden_Catchment.* Outline of the Study area and sub-catchments. Derived from gauge catchments in the NRFA
Eden_Urban.* Locations of urban areas as used in the 2011 census, from data.gov.uk
Eden_Catchment.* Outline of the Study area and sub-catchments. Derived from gauge catchments in the NRFA
Eden_River_Network.* River Network within the catchment. Taken from OS Open Rivers
Eden_DEM.tif Digital Elevation Model generated by aggregating OS Terrain 50 to 100m resolution
Eden_Precip.nc A time series of rainfall fields. Derived from the GPM IMERG final product by converting accrued rainfall in m

2.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

3 Initial GIS processing

3.1 Required R Packages

R has multiple packages for the processing and analysis spatial data. A good overview is given in the Spatial Task View on CRAN. For this part of the training course we use a set of mature package terra which is a dependency of dynatop so should already be installed. Similar results could be generated with other packages.

Attached the terra package to the R environment so that its functions are available.

library(terra)
library(dynatop)

and remember to move to the eden_data directory e.g.

setwd("../eden_data")

3.2 Data Required by dynatopGIS

To build a Dynamic TOPMODEL the dynatopGIS function requires the following

  • A vector polygon of the catchment outline, possibly delineating sub-catchments
  • A raster DEM covering the catchment
  • A vector line representation of the river network

Optionally the user can also provide

  • A raster map delineating spatial rainfall patterns
  • A raster map delineating spatial PET patterns
  • Raster maps containing further data to help classify the HRUs

All raster data must in the same resolution and projection as the dem. Vector data is automatically projected to match the dem.

dynatopGIS only works with projected raster data with square cells and assumes all raster values are numeric! The extent to the dem should be larger then the catchment outline so that no edge cells fall withing the catchment.

3.3 Loading Data

There are four types of spatial data in the example; raster, polygons, lines and points; in three different formats; Geotiff, shapefile and csv.

To start use functions in the raster package to load both the raster data and that contained in the shapefiles:

dem <- rast(file.path(".","unprocessed","Eden_DEM.tif")) # load the dem as a raster layer
eden <- vect(file.path(".","unprocessed","Eden_Catchment.shp")) # load the outline of the sub-catchments from the shapefile
channel <- vect(file.path(".","unprocessed","Eden_River_Network.shp")) # load the river network from the shapefile
urban <- vect(file.path(".","unprocessed","Eden_Urban.shp")) # load the urban area map from the shapefile

Typing the variable name at the R prompt will display a summary of the data loaded, for example

dem

Compare the dem, eden, channel and urban variable to check that:

  • there is a common projection
  • the cells of DEM are square

Reading the gauge locations from the csv file is a two stage process

 # read in the csv file to a date.frame
csvGauges <- read.csv(file.path(".","unprocessed","Eden_Gauge_Sites.csv"))
# convert the data frame to a SpatVect object - take projection from the DEM
gauges <- vect(csvGauges,geom=c("x","y"),crs=crs(dem))

Next plot the data using the default methods in the raster package

plot(dem) # underlying coloured image with scale
plot(urban,add=TRUE,col="grey") # grey polygons
plot(channel,add=TRUE,col="blue") # channels as blue lines
plot(eden,add=TRUE) # outlines of sub-catchments in black
plot(gauges,add=TRUE,col="red") # gauges as red dot

3.4 Rasterising Layers

The basis of the landscape discretisation is the raster catchment map. In this example, as in many situation the catchment outlines are vectors shapes. These are rasterised using the grid provided by the raster DEM. In the catchment map each sub-catchment requires a unique numeric value which is provided by the id field.

## show data frame with id values
head(eden)
## rasterise
eden_raster <- rasterize(eden, dem, field="id")
## plot the resulting raster
plot(eden_raster)

To ensure the edges are NA values the raster is extended, along with that of the DEM

eden_raster <- extend(eden_raster,1)
dem <- extend(dem,1)

Since the urban areas may be used in the classification they are also rasterised. Suitable numeric values are taken from the unique identifiers in the data.

In the case of the urban data there are two possible identifiers objectid and bua_id:

## show data frame with id values
head(urban)
## rasterise using objectid as the raster value
urban_raster <- rasterize(urban, eden_raster, field="objectid")
## plot the resulting raster
plot(urban_raster)

The raster fields created (eden_raster, urban_raster and extended dem) exist only in memory (or temporary files). Since they will be needed later we will save them:

writeRaster(eden_raster,file.path(".","processed","eden.tif"),overwrite=TRUE)
writeRaster(urban_raster,file.path(".","processed","urban.tif"),overwrite=TRUE)
writeRaster(dem,file.path(".","processed","dem.tif"),overwrite=TRUE)

3.5 Channel Network

The channel network in the dynamic TOPMODEL generated by dynatopGIS is a series of connected HRUs each representing a single reach length. In the processing done by dynatopGIS each reach is represented by a single spatial polygon. The data for each reach must include the following fields:

  • name - The unique identifier/name of the reach
  • length - Length of the reach [m]
  • startNode - An identifier for the head of the reach
  • endNode - An identifier for the foot of the reach
  • width - Channel width [m]
  • slope - Channel slope [m/m]

The startNode and endNode values are used to define the channel connectivity, while the physical properties are used in the dynatop numerical solutions.

Looking at the channel data

channel

we see that it already contains

  • an identifier for the foot of the reach called endNode
  • an identifier for the head of the reach called startNode
  • the reach length
  • a unique identifier called identifier

but is composed to spatial lines rather then polygons.

The dynatopGIS package provides a function for conveniently converting line objects to the correct format.

conv_channel <- convert_channel(
    channel,
    property_names = c(name = "identifier", length = "length", startNode = "startNode",
                       endNode= "endNode"),
    default_width = 2,
    default_slope = 0.001
)

The two warning messages reveal that, since they are not specified, default width and slope values are added. The width property is used to buffer the channel resulting in spatial polygons.

conv_channel

The converted channel data is saved for later use

writeVector(conv_channel,file.path(".","processed","channel.shp"),overwrite=TRUE)

The use of polygons to represent channel sections allows the shape of water features such as lakes to be captured

3.5.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

4 Observed series

4.1 Expected Data

The data passed to dynatop when simulating the model is expected to be a xts object (a special kind numeric matrix whose row names give the time) with a constant time step \(\Delta t\). Each column of data should have a unique name which is used to specify the time series when building the model.

Two variables are used as inputs, precipitation and potential evapotranspiration (PET). It is expected that the precipitation and potential evapotranspiration inputs series are given in \(m\) accrued over the proceeding time step. So the data value given at time \(t+\Delta t\) is the total accrued in the interval between time \(t\) and time \(t+\Delta t\).

Since columns in the xts object passed to dynatop which are not used in the model are ignored, so there is no danger is combining observed discharge records in the same object as the input data.

4.2 Data Files

The example contains two data files Eden_Flow.csv and Eden_Precip.nc. The first is a comma seperated text file of flow data, the later a netcdf file containing a time series of spatial rainfall fields.

4.3 Flow Data

Looking at the start of Eden_Flow.csv

cat(readLines(file.path(".","unprocessed","Eden_Flow.csv"),5),sep="\n")

we see that the Index column gives the time in an unabiguous format recognised by xts. The second column the observed discharges at the Sheepmount gauge as numeric values. We read this in using function avaialbe through the xts package

flow <- as.xts(read.zoo(file.path(".","unprocessed","Eden_Flow.csv"),header=TRUE,sep=",",drop=FALSE))
head(flow)

4.4 Rainfall Data

To explore the rainfall data open the file using the terra package

## read in as a multilayered object
precip_brk <- rast(file.path(".","unprocessed","Eden_Precip.nc"))
## show summary
precip_brk
## show start of the z values which are the time stamps as character strings
head(time(precip_brk))

From the properties of precip_brk we see that:

  • there are 1392 raster layers covering February 2020
  • the projection and resolution of the rainfall grids are different to the DEM
  • the time step is different to that of the flow data

To address the projection and resolution differences two approaches could be taken

  • reproject the rainfall data to match the DEM then extract a time series for each cell, so giving each DEM cell its own rainfall time series
  • assign each cell in the DEM to one of the rainfall data raster cells and extract the exisiting rainfall data to the give time series.

Since the DEM has a higher spatial resolution then the rainfall the second approach is more computationally efficent and is recommended.

To impliment the approach start by creating a raster layer of the same projection and extent as the rainfall feilds but with unique values

rid <- rast(precip_brk[[1]])
rid <- setValues(rid,1:ncell(rid))
rid

Then extract the precipitation values to give an xts object where each column is the time series for one value in rid.

precip <- values(precip_brk) # matrix where each row is one cell
## add row names based on unique values in rid to give time series names
rownames(precip) <-  paste0("precip_",values(rid))
## get times as characters and convert to R internal time representation (POSIXct)
precip_datetime <- time(precip_brk)
## create the xts object
precip <- xts(t(precip),order.by=precip_datetime)

To fix the difference in time step compared with the flow data we resample the precipitation data to a 15 minute timestep. The resample_xts function of dynatop provides a simple method for doing this

## call the resample_precip function
resampled_precip <- resample_xts(precip,900)
## look at what it does (assigned rainfall equally between new steps in same time period)
head(precip[,"precip_48"]) # original data
head(resampled_precip[,"precip_48"]) # resampled data - see how they sum to the original values

To assign the DEM cells to rainfall series reproject the raster layer of unique rainfall series numbers (rid) to match the projection and resolution of the DEM using nearest neighbours, then save this for later use.

## load the dem
eden <- rast(file.path(".","processed","eden.tif"))
## call the resample_precip function
reproj_rid <- project(rid,eden,method="near")
## A plot for visualisation
plot(reproj_rid)
## save
writeRaster(reproj_rid,file.path(".","processed","precip_id.tif"),overwrite=TRUE)

4.5 Potential Evapotranspriation

If available spatial Potential Evapotranspiration time series can be handled in exactly the same way as the rainfall. However if there is only a single series generatation of the map of id values can be avoided. A simple sinusoidal PET series can be generated using the evap_est function in dynatop.

## use a maximum of 3 mm/day - output in meters
pet <- evap_est(index(flow), eMin = 0, eMax = 0.003)
names(pet) <- "pet" # ensure series has a name

4.6 Combining Series

The xts package has a convient functions for merging xts objects by the time. We will save this as an R object for later use

obs <- merge(flow,resampled_precip,pet,all=FALSE) ## merge by time stamp
saveRDS(obs,file.path(".","processed","obs.rds"))

WARNING dynatop does not allow for missing data. Any missing values must be replaced during the pre-processing

5 Model Building

The dynatopGIS packages works through a number of steps to generate a dynamic TOPMODEL object suitable for use in with the dynatop package. Each step generates one or more layers which are saved as raster or shape files into the projects directory (which is not necessarily the R working directory).

Initialise the analysis by creating a new dynatopGIS object specifying the location of the project directory, which will be created if it doesn’t exist.

ctch <- dynatopGIS$new(file.path(".","dynaGIS"))

If existing files are present in the directory this step will try to reload them as a dynatopGIS project.

5.0.1 Adding catchment data

The basis of the analysis is a raster catchment map. Onto this a Digital Elevation Model (DEM) of the catchment and a vector representation of the river network with attributes will be added. Currently these can be in any format supported by the terrapackage.

However, within the calculations used for sink filling, flow routing and topographic index calculations the 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.

The catchment map needs to be loaded first

ctch$add_catchment( file.path(".","processed","eden.tif") )

Either the DEM or channel files can be added to the project next. In this case add the DEM processed during the initial GIS analysis:

ctch$add_dem(file.path(".","processed","dem.tif"))

As part of loading the DEM any blocks of NA values within the DEM that fall within the catchment are filled to appear as sinks. These will be filled in a later step.

Adding the channel data is achieved with

ctch$add_channel(file.path(".","processed","channel.shp"))

5.1 Getting and plotting catchment information

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

ctch$get_layer()

Note that channel_vect layer is the vector representation of the channel, while channel layer is a rasterised version.

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

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

or returned to the R workspace, for example

ctch$get_layer("dem")

5.2 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. In calling the sink_fill method a flow direction algorithm is specified and the resulting flow paths recorded. If sub-catchments are present in the catchment map then only flow paths within the sub-catchment are considered.

Since the algorithm 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 the Eden the DEM is already well conditions and the algorithm only alters a small areas of the catchment.

ctch$sink_fill()

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

5.3 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 1, those cells or reaches draining only into band 1 are in band 2 and so forth. Banding is achieved by the following call

ctch$compute_band()
ctch$plot_layer("band")

5.4 Computing Properties

Two sets of properties are required for Dynamic TOPMODEL. The first set is those required within the evaluation of the Hill slope HRU; gradient and contour length. The second set are those used for dividing the catchment up into different response classes. Traditionally the summary used for the separation of the classes is the topographic index, which is the natural logarithm of the up slope area divided by gradient.

These properties are computed using the formulae in Quinn et al. 1991.

The upstream area is computed by routing down slope with the fraction of the area being routed to the next downstream pixel being proportional to the gradient times the contour length.

The local value of the gradient is computed using the average of a subset of between pixel gradients. For a normal ‘hill slope’ cell these are the gradients to down-slope pixels weighted by contour length.

These properties are computed in an algorithm that passes over the data once in descending height. It is called as follows

ctch$compute_properties()

The plot of the topographic index shows a pattern of increasing values closer to the river channels

##plot of topographic index (log(a/tan b))
ctch$plot_layer('atb')

Although not used in ordering the HRUs dynatopGIS also provides the ability to compute flow distances for the hill slope cells. The calculation of three distances is supported

  • Shortest flow length - the shortest length based on the pixel flow paths to a channel
  • Dominant flow length - the distance to a channel moving in the dominant (largest fraction) flow direction from any grid cell
  • Expected flow length - the distance to the channel based on a weighted average of the down-slope flow lengths. Weights are given by the fraction of flow in each direction.

The computation, in this example for the shortest flow length, is initiated with

ctch$compute_flow_lengths(flow_routing="shortest")

The additional layers can be examined as expected

ctch$get_layer()
ctch$plot_layer("shortest_flow_length")

5.5 Adding additional Layers

Properties for use in the classification of HRUs or for determining input series may be added as additional raster GIS layers.

In the initial GIS analysis a layer for urban area was developed, while processing the input data produced a map relating to rainfall input. These can now be added to the dynatopGIS project by specifying the file location and name to call the layer.

ctch$add_layer(file.path(".","processed","urban.tif"),"urban")
ctch$add_layer(file.path(".","processed","precip_id.tif"),"precip_id")

Calling the get_layer method shows these are added to the project

ctch$get_layer()

5.6 Classifying into Hydrological Response Units

Methods are provided for the classification of the catchment. Multiple classifications can be produced but that used to generate the dynatop model (see following sections) must include the “band” layer.

By definition each channel length is treated as a single class with its own Channel HRU. Two ways are provided for dividing the hill-slope up into classes. The first way is cutting where a landscape property is divided up into classes. The second is combining where unique existing class layers are combined by either by pairing (where unique combinations of classes form a new class) or burning (which enforces classes onto distinct areas).

5.6.1 Cutting

To split a catchment into HRUs using cuts the breaks between the classes need to be specified. The breaks can be specified as either:

  • A single value: defines the number of splits which are automatically selected
  • A vector of values: treated as the breaks between classes

To demonstrate the use of cuts alone dividing the topographic index into 20 classes. Specifying a single value gives class breaks that are equally spaced across the range of the variable.

## simple equally spaced cuts
ctch$classify("atb_20","atb",20)

It is perhaps better to consider 20 classes with approximately equal land area. it this case the breaks can be computed using the quantiles of the variable

## use the quantile function of the raster package to compute breaks
## note we need to include the min and max
brk <- global(ctch$get_layer("atb"), fun=quantile,probs = seq(0, 1, length = 21),na.rm=T)
ctch$classify("atb_20_equal","atb",brk)

In the meta data the cuts and burns used to generate the layer are recorded. These can be retrieved with

head(ctch$get_method("atb_20"))
## returns the list of break points for the classes

5.6.2 Combining

The sub-catchments could be added alongside the banding to add a further spatial separation.

ctch$combine_classes("atb_20_subcatch",c("atb_20","catchment"))

If this is to be used to generate a model the band layer would have to be included.

ctch$combine_classes("atb_20_subcatch_band",c("atb_20_subcatch","band"))

A more complex combination would also include the urban area. Layers burnt into the classification are added after the cuts are applied and treated as directly defined HRUs. In the example it may be desirable to burn in the Urban areas since these will have a different hydrological behaviour.

Since the layers burnt in are applied with no alteration to their values these may clash with those already defined by the cuts. One method that nearly always avoid this is to burn in classes with negative values. Let us then adapt the urban layer to have negative values.

tmp <- -ctch$get_layer("urban")
writeRaster(tmp,file.path(".","processed","neg_extended_eden.tif"),overwrite=TRUE) ## write out
ctch$add_layer(file.path(".","processed","neg_extended_eden.tif"),"neg_urban")

then burn this into the existing classification based on the topographic index and sub-catchment

ctch$combine_classes("atb_20_subcatch_band_burn_neg_urban","atb_20_subcatch_band",burns="neg_urban")

In the meta data the classes of the original class values in the combination are stored These can be retrieved with

tmp <- ctch$get_method("atb_20_subcatch")
tmp$type
head(tmp$groups)
## returns the list of break points for the cuts and layers burnt in

5.7 Complexity

The number of classes resulting from the classifications reflects the complexity of the model. The number of classes can be computed from the maps and is shown in the following table. Since these classes may occur in any of the bands the resulting number of HRUs, which is also shown is much higher.

5.8 Generating a dynamic TOPMODEL

A Dynamic TOPMODEL suitable for use with the dynatop package can be generated using the create_model method. This uses an existing classification to generate the model.

For example, in the case of the division of by topographic index into 21 classes and the bands directly the resulting model can be generated by

ctch$create_model(file.path("./dyna","atb_20_model"), # name of new model
                  "atb_20_band", # classification to base model on
                  sf_opt="cnst",
                  sz_opt="exp",
                  rain_layer = "precip_id", # layer of input precipitation series ID
                  rain_label = "precip_" # characters added to values in rain_layer to get series name
                  )

Looking at the files within the folder containing the dynatopGIS project

list.files(file.path("./dyna"),pattern="atb_20_model")

shows that an additional raster map of the HRUs has been created in atb_20_model.tif along with a file atb_20_model.rds containing a model suitable for dynatop

5.8.1 Looking at the dynatopGIS model output

The model output by dynatopGIS should be adequate for dyantop to run, and is more fully documented here. Looking at the output in more detail shows atb_20_model.rds contains three variables.

## load the model
mdl <- readRDS( file.path(".","dyna","atb_20_model.rds")) ## read the model in
names(mdl)

The map variable contains the path of the HRU map.

mdl$map

The hru variable contains a description of the HRUs. This is given as an list, with each list element representing the properties of one HRU.

names(mdl$hru[[8]])

A summary of properties is given below, detailed descriptions are found in the dynatop vignettes

Name Class Description
id integer The unique ID of the HRU
states numeric Named numeric vector of states
properties numeric Named numeric vector of properties
sf list Description of the surface zone solution
rz list Description of the surface zone solution
uz list Description of the unsaturated zone solution
sz list Description of the saturated solution
sf_flow_direction list Description of the surface zone outflows to other HRUs
sz_flow_direction list Description of the saturated zone outflows to other HRUs
initialisation numeric Named numeric vector of initialisation parameters
precip list Description of the precipitation input to the HRU
pet list Description of the PET input to the HRU
class list Classification information from dynatopGIS

The output_flux variable contains a table in the format used in a call to a dynatop simulation to set the time series output.to set the

mdl$output_flux

By default dynatopGIS populates this to return the lateral outflows of the surface zones (flux column set to q_sf) of the catchment outlets (HRU number in the id column) giving then representative names.

This can be replaced by the locations of the river gauges found in the initial GIS processing.

## To replace this read in the gauge locations generated mealier
gauges <- vect(file.path(".","processed","gauges.shp"))
## create a vector of the uid's of the HRUs
hid <- sapply(mdl$hru,function(h){h$id})
## create a vector of the channel names with as storaged in the class variabel
cname <- sapply(mdl$hru,function(h){h$class$name})
idx <- match(gauges$chn_identi,cname) ## index of HRU's corresponding to gauges
## create as a data.frame
gauge_output <- data.frame(name = gauges$Site.Name,
                           id = as.integer(hid[idx]),
                           flux = "q_sf",
                           scale = 1)
saveRDS(gauge_output,file.path(".","dyna","gauge_output_defn.rds")) ## save the output for later use
gauge_output

5.8.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1

6 Model Simulation

To simulate a Dynamic TOPMODEL of the Eden catchment, visualise the output showing how the show how the parameters may be altered.

This exercise builds upon the early model creation, using the outputs generated during that exercise.

6.1 Locating Gauges on the Channel Network

Before locating the gauges on the network note that in the map of the loaded data there is a single gauge outside to the catchment area. As simple example of manipulating spatial data in R is to remove this by cropping the gauges to the catchment area. This is done by

gauges <- crop(gauges,eden)

To determine the location of the gauges on the river network we identify for each gauge the closest reach. The identifier of this reach, along with the distance to the gauge data is then added to the gauge data.

dst <- distance(gauges,conv_channel) ## distance between gauges(row) and channels(column)
idx <- apply(dst,1,which.min) ## index of closest channel for each gauge

gauges$chn_identifier <- conv_channel$name[ apply(dst,1,which.min) ]
gauges$chn_distance <- apply(dst,1,min) ## distance to the channel

To check the results by plotting:

selected <- conv_channel[conv_channel$name %in% gauges$chn_identifier,]
plot(eden) # outlines of sub-catchments in black
plot(conv_channel,add=TRUE,col="blue",border="blue") # channel network in blue
plot(selected,add=TRUE,col="orange",border="orange") # selected reachs in orange
plot(gauges,add=TRUE,col="red",pch=21) # gauges as red filled circles

The cropped list of gauges with channel identifier can be saved to a shapefile:

writeVector(gauges,file.path(".","processed","gauges.shp"),overwrite=TRUE)

The quality of results of any method used for any automatic method of locating the gauges depends upon the accuracy of the GIS data.

In the numeric solution used in dyantop it is assumed that a gauge is at the foot of a reach. This should be allowed from when conceptualising the channel reaches.

6.2 A Simulation

The minimum data need to perform a simulation with dynatop is a model, such as that output by dynatopGIS in the previous section and observed data in the format introduced earlier.

Load the previously created data

mdl <- readRDS(file.path(".","dyna","atb_20_model.rds"))
output_defn <- readRDS(file.path(".","dyna","gauge_output_defn.rds"))
obs <- readRDS(file.path(".","processed","obs.rds"))

A dynatop object is initialised by loading the model

hru <- mdl$hru
ctch_mdl <- dynatop$new(hru)

This step performs some basic checks on the model for consistency. The data can then be added; again simple checks are performed.

ctch_mdl$add_data(obs)

The initial states of the HRUs can be specified in the model object. If, as in the case of this example, the states are not specified then any attempt to perform a simulation will fail.

## expect a failure
ctch_mdl$sim(output_defn)

The states need to be initialised using the initialise method. This presumes:

  • The surface is dry
  • The root zone is partially full (the s_rz0 parameter)
  • There is a specified (see below) recharge rate representing the flow from the unsaturated to saturated zone
  • the saturated zone is at a steady state
  • the unsaturated zone storage ensures that the recharge rate from the unsaturated zone to the saturated zone is as close to the maximum as possible

The maximum recharge rate and fraction of root zone storage full for each HRU forms part of its definition.

hru[[1000]]$initialisation
## s_rz_0 - fraction of root zone full
## r_us_sz_0 - recharge unsaturated to saturated zone in m/s and

A simple way of getting an initial estimate of the recharge is to divide the observed discharge at the outlet of the catchment for the first time step by the catchment area. This works best if the first time step has no rainfall and is close to base-flow.

## determining the initial recharge rate ad setting it in the model
r0 <- as.numeric(obs$Sheepmount_obs[1]) / sum( sapply( hru, function(x){x$properties["area"]} ) )
hru <- lapply(hru,function(h){h$initialisation["r_uz_sz_0"] <- r0;h})

In the following we initialise the states and then perform the simulation, using the chaining of commands.

## the map input is only needed for visualisation
ctch_mdl <- dynatop$new(hru,map=mdl$map)$add_data(obs)$initialise()$sim(output_defn)

6.3 Mass balance

It is possible to output the mass balance for each time step of the simulation using the get_mass_errors method. The returns an xts object giving the volumes in the hill-slope states at the start and end of the time step along with the other fluxes as volumes. This can easily be used to plot the errors as shown below.

mb <- ctch_mdl$get_mass_errors()
head(mb)
plot( mb$err , main="Mass balance error",ylab="m^3")

6.4 Visualisations

After a simulation it is possible to get or plot output. Two types of output are available:

  • the time series variable specified in the output_defn used in the simulation
  • the states of the system at end of the simulation

Both have methods for plotting the data or returning in as a variable as shown in the following

6.4.1 Gauge flows

## For simulated flows at the gauge sites
sim_gauge <- ctch_mdl$get_output() ## extract the flow as an xts object
head(sim_gauge)

## Plot the simulated flows at the gauges
ctch_mdl$plot_output()

The observed gauge data is readily added and plotted

plot( merge(obs$Sheepmount_obs,sim_gauge),main="Simulated and observed flows",
     legend.loc='topright')

6.4.2 Current States

The current states of the model can be extracted and plotted.

## extract the current states as a data.frame
current_state <- ctch_mdl$get_states()
## columns correspond to the Hill-slope HRU ID, and the storage's in each of the four zones
head(current_state)

## a single state can be plotted at one time
ctch_mdl$plot_state("s_sz")

6.5 Changing Parameters

The parameters are stored within the HRU descriptions; particularly the sf, rz, uz and sz elements. Each of these contains a type variable and corresponding named vector of parameters. The various types and parameters for each zone are described in the (dynatop vignettes)[https://waternumbers.github.io/dynatop/articles/HRU.html#surface-zone-representations].

The most straight forward way for altering the parameters is two write a function that acts on each HRU, altering the parameters with reference to the properties of the HRU. The following presents a simple function of setting a faster surface velocity and kinematic solution for the river channels.

fpar <- function(h){ ## h is the hru
    ## test if a channel
    if( !is.na(h$class$startNode) ){
        ## change type
        h$sf$type <- "kin"
        h$sf$param["c_sf"] <- 1
    }
    h
}

6.6 Channel Inflow

To use dyantop to drive a hydraulic model the inflow from the subsurface and hillslope HRUs to the surface zone of each channel HRU may need to be computed. The first term is (directly available )[https://waternumbers.github.io/dynatop/articles/The_Model_Object.html#output_flux], though will the wrong sign and as a volume. Start by working out the id of each channel HRU

chn_id <- do.call(c, lapply( hru, function(h){ if(!is.na(h$class$startNode)){ h$id } }))

Next start creating an output definition for the flux from root zone

out_hyd <-  data.frame( name = paste0("hyd_",chn_id),
                       id = chn_id,
                       flux = "v_sf_rz", scale =-1/900) # scale is negative since flux usually downwards, 1/900 is 1/timestep

Next we will construct the inflow from HRUs representing hillslopes and add these to the data frame above.

## function to compute the relevent part so the output definition for each potential source HRU
fout <- function(h){
    ## test if it is a channel HRU - in which case return NULL
    if( !is.na(h$class$startNode) ){ return(NULL) }
        
    ## index of flow directiuons going to other channels
    idx <- h$sf_flow_direction$id %in% chn_id
    if(!any(idx)){ return(NULL) } ## quit if no links!
    ## create partial output defn
    data.frame(
        name = paste0("hyd_",h$sf_flow_direction$id[idx]), ## which inflow does it alter
        id = h$id, ## flux comes from this HRU
        flux = "q_sf", # surface flux out of this HRU is input to the next
        scale = h$sf_flow_direction$fraction[idx] # scale is fraction of flow (since taking it away)
    )
}
## row bind to give a data.frame
tmp <- do.call( rbind, lapply( hru, fout) )

## row bind onto the intial data frame
out_hyd <- rbind( out_hyd, tmp )

## sort this, not needed but easier to look at
out_hyd <- out_hyd[ order(out_hyd$name), ]

In the resulting output definition multiple fluxes contribute to the same output series, combined by a summation weighted by the scale value. Simulation and visualisation work as before.

ctch_mdl <- dynatop$new(hru)$add_data(obs)$initialise()$sim(out_hyd)
plot( ctch_mdl$get_output()[,1:10])

6.7 Extracting Intermediate States

Unless requested dynatop keeps only the final state during a simulation. To request that states are other time steps are kept a vector of times must be passed in the call to the sim method.

## Let us repeat the same simulation
## keeping the states for every other time step
keep_times <- index(obs)[seq(2,nrow(obs),by=96)] ## an vector of time steps
ctch_mdl$initialise()$sim(output_defn,keep_states=keep_times) ## passing the vector of time steps to the call

The record of the states can know be extracted. While there is an entry for every time step only those request will be populated

## extract the record of the states at each time step
state_rec <- ctch_mdl$get_states(record=TRUE) ## returns the whole record of states
head(names(state_rec)) ## the elements in state_rec are named after the time
head(sapply(state_rec,length)) ## those in keep_times are populated the remainder are empty

The visualisation of the state record is possible either as a time series

## Create a time series plot of the saturated zone deficit
s_sz_rec <- state_rec[ sapply(state_rec,length)>0 ]
s_sz_rec <- lapply(s_sz_rec, function(x){ setNames( sapply(x,function(z){z$states["s_sz"]}), sapply(x,function(z){z$id}) )} )
s_sz_rec <- do.call(rbind,s_sz_rec) ## bind to make matrix
s_sz_rec <- as.xts(s_sz_rec) ## convert to and xts object using the rownames as times
## for a the Hill-slope HRUs with IDs 2487, 3598 and 5201 the start of the time series are
head(s_sz_rec[,c("2487","3598","5201")])
## and a plot can be generated
plot(s_sz_rec[,c("2487","3598","5201")],main="Saturated zone deficit for select HRU",legend.loc="topright")

or as a map

## Create a map of the root zone storage at 2020-02-03 02:00:00
s_rz_val <- do.call(rbind, lapply( state_rec[["2020-02-03 00:30:00"]], function(x){ data.frame(id=x$id,s_rz = x$states["s_rz"]) }) )
## create a map of values by substituting s_rz for the id in the hru_id map
hru_id <- rast( file.path(".","dyna","atb_20_model.tif") ) ## open map from from file names stored in the model
s_rz_map<- subst(hru_id, from=s_rz_val$id, to = s_rz_val$s_rz) ## substitute values into the map of HRU id numbers
## and plot
plot(s_rz_map,main="Root zone storage at 2020-02-03 00:30:00")

6.8 Hot Starting and Saving

Hot starting; that is starting a simulation with the final states of a previous simulation can be achieved by simply not re-initialising the dynatop object before the next call to sim

## initialise and simulate the model
ctch_mdl$initialise()$sim(output_defn)
r1 <- ctch_mdl$get_output() ## get the output from the first run
## run again without initialising
## This will use the final states of the last simulation as the starting states
ctch_mdl$sim(output_defn)
r2 <- ctch_mdl$get_output() ## get the output from the second
plot(r1$Sheepmount)
lines(r2$Sheepmount,col="red")## different state values after the second run

If you want to save the model and current states at the end of a simulation do not save the dynatop object. Instead extract and save he model including the current states and save that.

new_mdl <- ctch_mdl$get_model() ## get the model structure from the dynatop object
new_mdl[[1]] ## as the original model but will final states included