The main purpose of the wateres package is calculation of various reservoir characteristics in order to evaluate its effectiveness. This tutorial describes the calculation steps and options for one single reservoir.

The package is focused on long-term reservoir balance in monthly time step, although some calculations for hourly and daily data are also supported.

Creating a wateres object

All input information about the reservoir is stored in an instance of the wateres class. To create this, a data frame with time series (currently date and average reservoir inflow in m³.s⁻¹) of input data is needed. Alternatively, a file containing the series can be used – a sample file of the Elven Rivendell is located within the tests of the package.

Together with the time series, reservoir potential storage and corresponding flooded area have to be specified:

library(wateres)
riv = as.wateres("../tests/testthat/rivendell.txt", 14.4e6, 754e3)

All package functions are implemented as methods of the wateres object (or data structures of outputs in case of plotting).

Setting additional input information

There are also optional time series and reservoir characteristics that affect the calculation.

For some water management purposes, water level outputs are required instead of volumes. Also, current flooded area of the reservoir could used when estimating its evaporation. Therefore, elevation-area-storage relationship can be set as an argument of the as.wateres function.

eas = data.frame(
    elevation = c(496, 499, 502, 505, 508, 511, 514, 517, 520, 523, 526, 529),
    area = 1e3 * c(0, 5, 58, 90, 133, 180, 253, 347, 424, 483, 538, 754),
    storage = 1e6 * c(0, 0.003, 0.161, 0.530, 1.085, 1.864, 2.943, 4.439, 6.362, 8.626, 11.175, 14.4))
riv = as.wateres("../tests/testthat/rivendell.txt", 14.4e6, 754e3, eas = eas)

Evaporation from the flooded area of the reservoir can be set directly as time series (of the length equal to the length of the inflow data or 12 monthly values; in mm) or it can be estimated by a method according to the Czech Technical Standard, where evaporation is a function of reservoir altitude.

All modification of the data are performed by a method of the wateres object -- the reservoir is passed as its first argument and a modified object is returned finally. Hence, the evaporation series for the Rivendell reservoir are set as follows:

riv = set_evaporation(riv, altitude = 529)

and can be accessed as riv$E.

Similarly, time series of water use from the reservoir are set. The water use values can be both positive and negative, meaning water release to the reservoir (added to the water balance in any case) or withdrawal from the reservoir (i.e. water demand of lower priority than yield), respectivelly. The set_wateruse function accepts a constant value, 12 monthly values or complete time series (in m³); the values are stored in riv$W.

riv = set_wateruse(riv, c(7, 14, 40, 62, 82, 96, 109, 102, 75, 48, 34, 13))

Finally, precipitation values (in mm) stored in riv$P are supported:

riv = set_precipitation(riv, c(55, 40, 44, 43, 81, 72, 85, 84, 52, 54, 48, 58))

Precipitation is supposed to affect the reservoir flooded area corresponding with the maximum storage.

Calculating reservoir water balance

To calculate time series of water balance variables of the reservoir, use the calc_series function:

resul = calc_series(riv, yield = 0.14, get_level = TRUE)

A required yield is the only argument to be set (as a fixed yield or a vector of values), however both maximum and initial storage can be also specified. Additionally, water levels can be calculated if elevation-area-storage relationship has been given and the get_level argument is set.

The output is returned as a wateres_series object, i.e. a data table with water balance variables. It can be easily visualized by its plot function using the ggplot2 package. Three plot types -- flows, storage and levels are supported:

p = plot(resul, riv, "flow", begin = 0, end = 500)
p = plot(resul, riv, "storage", begin = 0, end = 500)
p = plot(resul, riv, "level", begin = 0, end = 500)

Inspecting storage-reliability-yield relationship

To optimize reservoir storage with respect to the required yield and the related time-based reliability, use the sry function. This function is passed by values of two of these three characteristics (storage, yield and reliability; potential reservoir storage is considered as a default storage value) while the remaining one will be:

For instance, the reliability value for the yield required as above is found by:

sry(riv, yield = 0.14)

meaning that the yield value is ensured in each time step.

If a lower reliability is sufficient, a corresponding reservoir storage can be calculated as:

sry_resul = sry(riv, reliability = 0.95, yield = 0.14)
sry_resul

The calculated time series can be obtained by setting the get_series argument. There are also arguments controlling a method of reliability calculation^[Therefore, maximum reliability could be less than 1 and "max" value is also allowed as the reliability argument.] or bisection limits.

Calculating reservoir characteristics

The summary function is a more convenient and concise way how to obtain reservoir characteristics. It employs the sry function and requires the same arguments. Moreover,

summary(riv, reliability = c(0.9, 0.95), yield = 0.14)

Plot outputs

Time series of water balance components can be also analyzed in terms of their probability for particular months. The probability fields are supported for reservoir storage, level and yield. They are calculated and visualized by:

prob_field = prob_field(riv, probs = c(0.9, 0.95, 0.99), yield = 0.14, storage = sry_resul$storage)
p = plot(prob_field, "storage")
p = plot(prob_field, "yield")

There is also a traditional method to characterize the reservoir by using curves representing the reservoir effectiveness as relation between level of development (alpha, determining the yield value) and ratio of storage and volume of annual flow (beta). These curves can be displayed for a vector of reliability values:

ab = alpha_beta(riv, alphas = seq(0, 1, 0.1), reliability = c(0.9, 0.95, 0.99))
p = plot(ab)

Calculating catchment and reservoir deficit volumes

Estimation of reservoir effectiveness with respect to water demands in the corresponding catchment is a typical use case for the wateres package. Suppose that we need to cover water use requirements in the Bruinen catchment by the Rivendell reservoir. The time series of the Bruinen river discharge for the catchment outlet and water use within the catchment are available in a sample file (with a distinct change in water use; positive water use means that the river is supplied by water release originating e.g. in groundwater resources):

bru_data = read.table("../tests/testthat/bruinen.txt", header = TRUE, colClasses = c("Date", "numeric", "numeric"))
library(ggplot2)
ggplot(bru_data, aes(x = as.Date(DTM), y = USE)) + geom_line()

The catchment can be represented by a wateres reservoir whose storage is set to zero. There should be also a minimum residual flow set for the river:

bru = as.wateres(data.frame(DTM = bru_data$DTM, Q = bru_data$Q), storage = 0, area = 0)
bru_mrf = 1.514

To obtain water demand deficits, yield has to be set to the minimum residual flow increased by withdrawal requirements (negative water use) or decreased by water release (positive water use).^[Water use variable cannot be used in this case because withdrawal may occur only if there is a storage greater than zero.]

bru_resul = calc_series(bru, yield = bru_mrf - bru_data$USE, throw = TRUE)
p = plot(bru_resul, bru, "flow")

The remaining river flow cannot be visualized directly by the plot.wateres_series function, since it consists of the MRF (part of yield) and the difference between inflow and yield^[Specifically, the yield calculated with the throw argument set to TRUE, otherwise yield would equal to inflow because of no storage.]:

ggplot(
    reshape2::melt(
        data.frame(DTM = bru$DTM, Q_original = bru$Q, Q_with_use = bru_mrf + bru$Q - bru_resul$yield),
        id = "DTM"),
    aes(x = DTM, y = value, colour = variable)) + geom_line()

Afterwards, calculation for the reservoir can be performed. At first, a wateres object for the corresponding shorter time series is created -- the time series of the previous riv object can be shortened by the resize_input command which keeps the evaporation and precipitation values:

idcs = which(riv$DTM %in% bru$DTM)
riv_short = resize_input(riv, idcs[1], idcs[length(idcs)])

Then, set withdrawal (i.e. negative water use) equal to the catchment deficits and a minimum residual flow value for the reservoir profile.^[The use of this MRF produces a bias since minimum residual flow demand has been already included in the catchment deficits; however it is applied to make calculations for catchment and reservoir consistent.]

riv_short = set_wateruse(riv_short, -bru_resul$deficit)
riv_mrf = 0.033

Having catchment deficits defined as water use, reservoir yield is equal to the minimum residual flow:

riv_resul = calc_series(riv_short, yield = riv_mrf)
p = plot(riv_resul, riv_short, "flow")
p = plot(riv_resul, riv_short, "storage")

implying that the deficits would be reduced by the reservoir by 1 - sum(riv_resul$deficit) / sum(bru_resul$deficit), i.e. by r 100 * round(1 - sum(riv_resul$deficit) / sum(bru_resul$deficit), 2) %.

Using hourly data

Use of the wateres package is not limited to the long-term water balance estimation. Since the package supports also hourly and daily time steps and variable yield output, it can be also applied to e.g. flood wave transformation:

rivh = as.wateres("../tests/testthat/rivendell_1h.txt", 14.4e6, 754e3, eas = eas, time_step = "hour")
resul = calc_series(
    rivh, yield = c(rep(0.2, 8), rep(1, 16)), get_level = TRUE, initial_storage = 14e6)
p = plot(resul, rivh, "flow")
p = plot(resul, rivh, "storage")


tgmwri/wateres documentation built on Feb. 13, 2024, 10:25 p.m.