Get started with transfR
In transfR: Transfer of Hydrograph from Gauged to Ungauged Catchments

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This package aims to estimate discharge time series for ungauged catchments (non-instrumented catchments where discharge is not available) using hydrological observations from neighbouring gauged catchments (instrumented catchments where discharge is available). The hydrological modelling is based on a description of catchment geomorphology that can be assessed at any location. Inverting this model at gauged locations makes it possible to estimate net rainfall, which facilitates the transfer of observed discharge to ungauged locations.

1. Create a transfR object

An object of class transfR must first be created with the as_transfr() function. It is also used to gather all catchment attributes and intermediate results from the different steps. This object requires two user inputs:

a spatio-temporal object (a stars object) describing both discharge and the corresponding spatial support (outlet, centroid or catchment boundary); see the vignette Preparation of input data: creation of a stars object;
raster maps of hydraulic length (stars or matrix objects) for each catchment, describing the flow path length from each pixel to the outlet within the river network [@Cudennec2004; @Aouissi2013]; see the vignette Preparation of input data: geomorphological analysis with whitebox.

This package does not provide functions to create these two inputs. They must be prepared beforehand by the user. Several GIS software packages can extract them from a digital elevation model, including GRASS toolkits [@Jasiewicz2011], Whitebox GAT (see @Lindsay2016 or WhiteboxTools), TauDEM (D. Tarboton, Utah State University) and online services (@Squividant2015 for catchment delineation only). The vignettes mentioned above provide guidance on preparing the input data.

The 'Oudon' example dataset contains these two inputs, with hourly discharge observations for six sub-catchments of the Oudon River in France, their respective catchment boundaries and maps of their hydraulic lengths. All catchments are gauged; however, in this example, we use the first three as gauged catchments and the last three as ungauged catchments. To evaluate the methodology, the package can also perform a leave-one-out analysis by treating each gauged catchment in turn as ungauged, without creating the sim object below (set cv=TRUE when using mixr() in step 6).

library(transfR)
data(Oudon)
obs <- as_transfr(st = Oudon$obs[,,1:3], hl = Oudon$hl[1:3]) #gauged catchments
sim <- as_transfr(st = Oudon$obs[,,4:6], hl = Oudon$hl[4:6]) #catchments considered as ungauged

The package also offers simple plots for transfR objects.

plot(x = obs, i = 1, attribute = "Qobs", format = "%b %d")
plot(obs$hl[[1]], axes = T, main = "Hydraulic length of gauged catchment 1 [m]",
     downsample=1,col=hcl.colors(n=20,palette="Blues"))

2. Streamflow velocity

The streamflow velocity (uc) is the unique parameter of the transfer function that needs to be estimated. It allows assessing the travel time from each pixel of the map of hydraulic length (hl) to the outlet, and then to create the unit hydrograph (uh). Here we will use the function velocity() to estimate an average streamflow velocity from a regionalisation established over the Loire River [@deLavenne2016; see help of the function for details]. If the input of velocity() is a transfR object, the velocity will be computed for each catchment.

obs <- velocity(obs, method = "loire2016")
obs$uc
sim <- velocity(sim, method = "loire2016")
sim$uc

3. A geomorphology-based unit hydrograph

Assuming a description of the flow path length (hl) and a streamflow velocity (uc), the transfer function of the river network can be built based on unit hydrograph theory. This is done using the function uh(), which requires these two inputs. If the input of uh() is a transfR object, the unit hydrograph will be computed for each catchment.

obs <- uh(obs)
sim <- uh(sim)
plot(obs, i = 1, attribute = "uh")

4. Net rainfall estimated a priori

To solve the inversion, an a priori estimate of net rainfall must be provided. We estimate it from the specific discharge, shifted by a lag time using the lagtime() and rapriori() functions.

obs <- lagtime(obs)
obs <- rapriori(obs)

5. Net rainfall estimated by inversion

Using the results of the previous steps, the inversion can be computed to estimate the net rainfall time series for each gauged catchment that best reproduces the observed discharge according to the transfer function. This follows inversion theory [@Tarantola1982; @Menke1989; @Boudhraa2018; see the help for inversion() for more details]. If the input to inversion() is a transfR object, net rainfall time series are estimated for each gauged catchment sequentially by default, or across catchments in parallel when parallel = TRUE. Parallel processing can reduce computation time when the object contains several catchments with long time series.

obs <- inversion(obs)

6. Estimate net rainfall at ungauged locations

The net rainfall of an ungauged catchment is estimated by averaging the net rainfall of neighbouring gauged catchments. This average can be weighted by the inverse of the distance between each gauged catchment and the ungauged catchment. The distance between two catchments is the rescaled Ghosh distance, computed with hdist() as defined by @deLavenne2016. The mixr() function then uses this distance matrix to estimate the net rainfall time series at every ungauged location.

mdist <- hdist(x = obs, y = sim, method = "rghosh", parallel = TRUE, cores=2)
sim <- mixr(obs = obs, sim = sim, mdist = mdist)

7. Simulate discharge at ungauged locations

Discharge time series at ungauged locations can finally be simulated through a convolution between the unit hydrograph and the net rainfall time series of each catchment. This is done with the convolution() function. In this example, simulated and observed discharge can be compared because the ungauged locations were deliberately defined from gauged locations (as described in step 1). Note that the beginning and end of the simulation are removed because of the warmup and cooldown periods required for inversion.

sim <- convolution(sim)
plot(x = sim, i = 1, attribute = c("Qobs","Qsim"), 
     ylab = expression(paste("Discharge [",m^3/s,"]")),
     col = c("#a6bddb","#045a8d"), format = "%b %d")