Examples DWB model

Camila Garcia and David Zamora

2026-07-17

DWBmodelUN examples in a tropical basin

DWBCalculator function

DWBCalculator is the function responsible for incorporating the DWB model calculations into the DWBmodelUN package. In this connection, the following example corresponds to an ungauged modeling exercise developed for the Sogamoso River basin (SRB) by (Duque, 2018).

As a first step, the example starts with loading datasets generated for SRB that are required by DWB model and subsequently, the coordinates system of the input data are verified. The last step is to run DWBCalculator function in order to obtain a monthly streamflow simulated.

library(DWBmodelUN)
#> Cargando paquete requerido: raster
#> Cargando paquete requerido: sp
# Load P and PET databases
data(P_sogamoso, PET_sogamoso)
# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)
#> First data file is a data frame
#> Second data file is a data frame
#> Two data frames - Comparing coordinates and headers
#> Coordinates verified
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): First date/header does
#> not match - Please verify
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): Final date/header does
#> not match - Please verify
#> Extent verified
#> Resolution verified
#> Rows and columns verified
#> Number of layers verified
# Load geographic info of GRU and parameters per cell
data(GRU, param)
# Construction of parameter maps from values by GRU
GRU.maps <- buildGRUmaps(GRU, param)
alpha1_v <- GRU.maps$alpha1
alpha2_v <- GRU.maps$alpha2
smax_v <- GRU.maps$smax
d_v <- GRU.maps$d

# Establish the initial modeling conditions
init <- init_state(GRU.maps$smaxR)
g_v <- init$In_ground
s_v <- init$In_storage
rm(init)

# Load general characteristics of modeling
setup_data <- readSetup(Read = TRUE)
Dates <- seq(as.Date( gsub('[^0-9.]','',colnames(P_sogamoso)[3]), format = "%Y.%m.%d"), 
             as.Date(gsub('[^0-9.]','',tail(colnames(P_sogamoso),1)) , format = "%Y.%m.%d"), by = "month")
Start.sim <- which(Dates == setup_data[8,1]); End.sim <- which(Dates == setup_data[10,1])
# Sim.Period: the 1st two columns of the P and PET are the coordinates of the cells
Sim.Period <- c(Start.sim:End.sim)+2  

# Run DWB model
DWB.sogamoso <- DWBCalculator(P_sogamoso[ ,Sim.Period], 
                              PET_sogamoso[ ,Sim.Period],
                              g_v, s_v, alpha1_v, alpha2_v, smax_v, d_v)

Calibration DWB model with Dynamical Dimension Search algorithm

dds is the function that implements the calibration algorithm Dynamical Dimension Search by (Tolson and Shoemaker, 2007). In the example, the DWB model parameters are calibrated in the SRB with 10 Group Response Unit (GRU). The Nash-Sutcliffe efficiency (nse.cof) is created as an objective function between simulated and observed streamflow data, with the treatment of missing values and its coupled with the DWBCalculator function. Finally, the calibration is done with 2 iterations using the databases of the sample basin.

library(DWBmodelUN)
# Load P and PET databases
data(P_sogamoso, PET_sogamoso)

# Verify that the coordinates of the databases match
Coord_comparison(P_sogamoso, PET_sogamoso)
#> First data file is a data frame
#> Second data file is a data frame
#> Two data frames - Comparing coordinates and headers
#> Coordinates verified
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): First date/header does
#> not match - Please verify
#> Warning in Coord_comparison(P_sogamoso, PET_sogamoso): Final date/header does
#> not match - Please verify
#> Extent verified
#> Resolution verified
#> Rows and columns verified
#> Number of layers verified

# Load geographic info of GRU and basins where calibration will be performed
data(GRU,basins)
cellBasins <- cellBasins(GRU, basins)
#> Warning in cellBasins(GRU, basins): gruLoc and basins have different coordinate
#> reference systems. The basins are reprojected to the raster CRS during the
#> extraction, but verify projections before using the results.
#> Warning: [extract] transforming vector data to the CRS of the raster

# Establish the initial modeling conditions
GRU.maps <- buildGRUmaps(GRU, param)
init <- init_state(GRU.maps$smaxR)
g_v <- init$In_ground
s_v <- init$In_storage
rm(init)

# Load general characteristics of modeling
setup_data <- readSetup(Read = TRUE)
Dates <- seq(as.Date( gsub('[^0-9.]','',colnames(P_sogamoso)[3]), format = "%Y.%m.%d"), 
             as.Date(gsub('[^0-9.]','',tail(colnames(P_sogamoso),1)) , format = "%Y.%m.%d"), by = "month")

# For this calibration exercise, the last date of simulation is 
# the same as the final date of calibration
Start.sim <- which(Dates == setup_data[8,1])
End.sim <- which(Dates == setup_data[11,1])
# the first two columns of the P and PET are the coordinates of the cells
Sim.Period <- c(Start.sim:End.sim)+2 
Start.cal <- which(Dates == setup_data[9,1])
End.cal <- which(Dates == as.Date("2004-12-01"))
# the first two columns of the P and PET are the coordinates of the cells
Cal.Period <- c(Start.cal:End.cal)+2  

#Load observed runoff
data(EscSogObs)

# Function that runs the DWB model
NSE_Sogamoso_DWB <- function(parameters, P, PET, g_v,s_v, Sim.Period, EscObs, Cal.Period){
  
  parameters <- as.vector(parameters)
  # Transform the parameters to the format that the model needs
  nGRU <- max(terra::values(terra::rast(GRU)), na.rm = TRUE)
  param <- matrix(parameters, nrow = nGRU)
  
  # Construction of parameter maps from values by GRU
  GRU.maps <- buildGRUmaps(GRU, param)
  alpha1_v <- GRU.maps$alpha1
  alpha2_v <- GRU.maps$alpha2
  smax_v <- GRU.maps$smax
  d_v <- GRU.maps$d
  DWB.sogamoso <- DWBCalculator(P_sogamoso[ ,Sim.Period], PET_sogamoso[ ,Sim.Period],
                                g_v,s_v, alpha1_v, alpha2_v, smax_v,d_v, calibration = TRUE)
  Esc.Sogamoso <- varBasins(DWB.sogamoso$q_total, cellBasins$cellBasins)
  
  # model evaluation; in case of possible NA results in the simulation,
  # add a conditional assignment to a very high value
  sim <- Esc.Sogamoso$varAverage[Cal.Period - 2, ]
  # align the observations with the simulated basins (named by gauge code)
  obs <- EscSogObs[Cal.Period - 2, paste0("X", colnames(sim))]
  
  if (sum(!is.na(sim)) == prod(dim(sim))){
    numer <- apply((sim - obs)^2, 2, sum, na.rm = TRUE)
    demom <- apply((obs - apply(obs, 2, mean, na.rm = TRUE))^2, 2, sum, na.rm = TRUE)
    nse.cof <- 1 - numer / demom
  } else {
    nse.cof <- NA
  }
  
  Perf <- (-1)*nse.cof
  if(!is.na(mean(Perf))){ 
    Mean.Perf <- mean(Perf)
  } else {Mean.Perf <- 1e100}
  return(Mean.Perf)
}

# coupling with the DDS algorithm
xBounds.df <- data.frame(lower = rep(0, times = 40), upper = rep(c(1, 2000), times = c(30, 10)))
result <- dds(xBounds.df = xBounds.df, numIter=2, OBJFUN=NSE_Sogamoso_DWB,
              P = P_sogamoso, PET = PET_sogamoso, g_v = g_v, s_v = s_v, Sim.Period = Sim.Period, 
              EscObs = EscSogObs, Cal.Period = Cal.Period)

Interactive graphics to represent outputs from the DWB model

graphDWB is the function in charge of creating the graphical inputs and outputs of DWBmodelUN. The function has four types of graphs to be generated:

  1. The first one (tp = 1) corresponds to a hydrograph of any continuous variable, such as the precipitation entries.
  2. The second one (tp = 2) compares the simulated runoff with the observations.
  3. The third one (tp = 3) generates a combined chart with a rain hyetograph and the simulated and observed streamflow hydrographs.
  4. The fourth one (tp = 4) adds a hydrograph of evapotranspiration to the previous comparison.

The present example shows how to generate three of the four chart types, using the datasets included in DWBmodelUN for the SRB. To draw one of the chart types it is necessary to load the variables involved in each one, to adjust the input series to a time vector and assign them to a list element. Finally to run graphDWB.

First: Example with precipitation time series

library(DWBmodelUN)
library(dygraphs)
#> Warning: package 'dygraphs' was built under R version 4.6.1
data(P_sogamoso)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
var <- list("Precipitation" = P.est)
 graphDWB(var, tp = 1, main = "Precipitation Lat:7.0 Lon:-72.94")

Second: Example with precipitation and runoff time series

library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, simDWB.sogamoso, EscSogObs)
P.est <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[, 1], start = c(2001, 1), frequency = 12)
runoff.obs <- ts(EscSogObs[, 1], start = c(2001, 1), frequency = 12)
var <- list("Precipitation" = P.est,"Runoff.sim" = runoff.sim, "Runoff.obs" = runoff.obs)
graphDWB(var, tp = 3, main = "DWB results at Sogamoso Basin")

Third: Example with precipitation and runoff time series

library(DWBmodelUN)
library(dygraphs)
data(P_sogamoso, PET_sogamoso, simDWB.sogamoso)
P <- ts(c(t(P_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
PET <- ts(c(t(PET_sogamoso[1, -2:-1])), start = c(2001, 1), frequency = 12)
runoff.sim <- ts(simDWB.sogamoso[, 1], start = c(2001, 1), frequency = 12)
var <- list("P" = P,"PET" = PET, "Runoff.sim" = runoff.sim)
graphDWB(var, tp = 4, main = "General Comparison Sogamoso Basin")

References

Duque, Nicolás. 2018. “Estimación de Campos de Precipitación En Cuencas Hidrográficas Colombianas Con Escasez de Datos, Combinando Datos Teledetectados Y de Estaciones En Tierra, Utilizando Funciones de Kernel.” Master’s thesis, Universidad Nacional de Colombia - Sede Bogotá. https://repositorio.unal.edu.co/handle/unal/71663.

Tolson, Bryan A, and Christine A Shoemaker. 2007. “Dynamically dimensioned search algorithm for computationally efficient watershed model calibration” 43: 1–16.