Páramo.be - research in the páramo ecosystem

Home
Overview
Monitoring
Projects
CV
Publications
R code
Links

Running TOPMODEL from R

For an overview see ?topmodel or consult the html help page that should be installed with the package

Here an example is given of a complete topmodel analysis, including presentation of the results. The most crucial part is to get the data in the right format, which of course depends on the format they come in. But with some understanding of R, it should not be too difficult.

Import and preprocessing input data

All data should be in simple vector format. The section on data processing section may be helpful.

The delay function

The delay function is used to calculate the residence time in the stream network, and thus to disperse both overland flow and subsurface flow generated over time.
The model uses a 2 column matrix. The first column represents cumulative fractions of the catchments (that sum to 1) and the second column represents the average distance of the stream network from each fraction to the outlet of the catchment (in m). Parameter vr is used to convert the delay function to the time domain.

A realistic delay function for a 1 km2 catchment could be:

fractions <- c(0.1,0.5,1)
distance <- c(200,500,700)
delay <- cbind(fractions, distance)

If the gauging location is further downstream than the catchment outlet, say 1000m, then the delay function would look like:

fractions <- c(0,0.1,0.5,1)
distance <- c(1000,1200,1500,1700)
delay <- cbind(fractions, distance)

Running the model

First, an array of parameters should be constructed. Currently, the model takes the following parameters:

qs0 = Initial subsurface flow per unit area (m)
lnTe = Areal average of ln(T0), the log of the saturated transmissivity (m/h).
m = Model parameter
Sr0 = Initial root zone storage deficit (m)
Srmax = Maximum root zone storage deficit (m)
td = Unsaturated zone time delay per unit storage deficit. Depending on whether this value is positive or negative, a different formula is used (see Beven et al., 1995).
vch = channel flow outside of the catchment. This parameter is only used if the gauging location is further downstream than the catchment outlet (m/h).
vr = channel flow inside catchment (m/h).
k0 = Surface hydraulic conductivity. This parameter is used in the unsaturated zone time delay function if td is negative.
psi = Wetting front suction for the calculation of infiltration excess flow
dtheta = Water content change across the wetting front.
dt = timestep (h)

Initialising the parameters, and running the model:

qs0 <- 2.662e-05
lnTe <- -0.56
m <- 1.336e-02
Sr0 <- 1.29e-02
Srmax <- 3.5e-03
td <- 0.75
vch <- 1000
vr <- 1134
k0 <- 10
psi <- 0.1
dtheta<- 1.3
dt <- 0.25

params <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,psi,dtheta,dt)

topmodel(params, topidxclasses, delay, rain, ET0, Qobs = Qobs)

NOTE: In most cases it is desirable to have the model return the simulated discharge, which can then be evaluated with an objective function. However, this evaluation (which can be done in R) may be slow. To speed things up, topmodel() has a built in shortcut to calculate and return the Nash - Sutcliffe efficiency directly. This happens when an array with observed discharges is passed to the function, as in the example above (Qobs). If Qobs is NOT given, an array with simulated discharges is returned.

Postprocessing

R provides many functions for postprocessing and visualisation. For instance:

Qsim <- topmodel(params, topidxclasses, delay, rain, ET0 )
# invoked without the Qobs argument, the simulated flow is returned
plot(Qobs, type="l")
points(Qsim, type="l", col="red")

References

Beven, K. J., and Kirkby, M. J. A physically based variable contributing area model of basin hydrology. Hydrol. Sci. Bull. 24 (1979), 43-69.
Beven K, Lamb R, Quinn P, Romanowicz R, Freer J, 1995. TOPMODEL. In: Sing VP (Ed), Computer Models of Watershed Hydrology. Water Resources Publications, Colorado. pp. 627-668.