Chapter 4: Semi-distributed glacio-hydrological modeling

library(HBV.IANIGLA)

data("glacio_hydro_hbv")

str(glacio_hydro_hbv)
#> List of 5
#>  $ baisn:'data.frame':   15 obs. of  7 variables:
#>   ..$ elev_band      : chr [1:15] "eb_1" "eb_2" "eb_3" "eb_4" ...
#>   ..$ soil_area(km2) : num [1:15] 10 15 17 20 25 18 17 16 14 12 ...
#>   ..$ ice_area(km2)  : num [1:15] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ total_area(km2): num [1:15] 10 15 17 20 25 18 17 16 14 12 ...
#>   ..$ rel_soil       : num [1:15] 0.0526 0.0789 0.0894 0.1052 0.1314 ...
#>   ..$ rel_ice        : num [1:15] 0 0 0 0 0 0 0 0 0 0 ...
#>   ..$ h(masl)        : num [1:15] 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 ...
#>  $ tair : num [1:5310, 1:15] 24.4 26.9 27 24.5 22.1 ...
#>  $ prec : num [1:5310, 1:15] 0 0 0 0 0 0 0 0 0 0 ...
#>  $ pet  : num [1:5310, 1:15] 0 0.999 0.997 0.993 0.988 ...
#>  $ qout : num [1:5310, 1:3] 0.0021 0.0091 0.0218 0.0438 0.0617 0.0711 0.0774 0.0809 0.0858 0.0936 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr [1:3] "total" "soil" "glacier"
#>  - attr(*, "comment")= chr "Semi-distributed glacio-hydrological HBV model input data"

## brief arguments description
  # basin: data frame with the same structure of the data("glacio_hydro_hbv) (colnames included).
  # tair: numeric matrix with air temperature inputs. 
  # precip: numeric matrix with precipitation inputs. 
  # pet: numeric matrix with potential eavapotranspiration inputs. 
  # param_snow: numeric vector with snow module parameters.
  # param_ice: numeric vector with glacier parameters.
  # param_soil: numeric vector with soil moisture parameters.
  # param_route: numeric vector with the routing parameters.
  # param_route_ice: numeric vector with the glacier routing parameters.
  # param_tf: numeric vector with the transfer function parameter.
  # init_snow: numeric value with initial snow water equivalent. Default value being 20 mm.
  # init_ice: numeric value with initial snow water equivalent of the glaciers. Default value
  # being 20 mm.
  # init_soil: numeric value with initial soil moisture content. Default value being 0 mm.
  # init_route: numeric vector with bucket water initial values. Default values are 0 mm.
  # init_route_ice: numeric value with glacier bucket initial value. Default values are 0 mm.
## output
  # simulated streamflow series.
glacio_hydrological_hbv <- function(basin,
                                    tair,
                                    precip,
                                    pet,
                                    param_snow,
                                    param_ice,
                                    param_soil,
                                    param_route,
                                    param_route_ice,
                                    param_tf,
                                    init_snow = 20,
                                    init_ice = 20,
                                    init_soil = 0,
                                    init_route = c(0, 0, 0),
                                    init_route_ice = 0
                                    ){
  n_it <- nrow(basin)

  # create output lists
  snow_module   <- list()
  ice_module    <- list()
  soil_module   <- list()
  route_module  <- list()
  route_ice_mod <- list()
  tf_module     <- list()

  # snow and soil module in every elevation band
  for(i in 1:n_it){
    snow_module[[ i ]] <-
      SnowGlacier_HBV(model = 1, inputData = cbind(tair[ , i], precip[ , i]),
                      initCond =  c(init_snow, 2), param = param_snow)

    ice_module[[ i ]] <-
      SnowGlacier_HBV(model = 1, inputData = cbind(tair[ , i], precip[ , i]),
                      initCond =  c(init_ice, 1, basin[i, 'rel_ice']), param = param_ice)

    soil_module[[ i ]] <-
      Soil_HBV(model = 1, inputData = cbind(snow_module[[i]][ , 5] , pet[ , i]),
               initCond = c(init_soil, basin[i, 'rel_soil']), param = param_soil )

  } # end for

  # get total soil discharge
  soil_disch <- lapply(X = 1:n_it, FUN = function(x){
    out <- soil_module[[x]][ , 1]
  })
  soil_disch <- Reduce(f = `+`, x = soil_disch)

  # get swe and total ice melt for all glacier area
  ice_disch <- lapply(X = 1:n_it, FUN = function(x){
    out <- ice_module[[x]][ , 9]
  })
  ice_disch <- Reduce(f = `+`, x = ice_disch)

  ice_swe   <- lapply(X = 1:n_it, FUN = function(x){
    out <- ice_module[[x]][ , 3] *  (basin[x, 'rel_ice'] / sum(basin[ , 'rel_ice']) )
  })
  ice_swe <- Reduce(f = `+`, x = ice_swe)

  # route module
  route_module <- Routing_HBV(model = 1, lake = F, inputData = as.matrix(soil_disch),
                              initCond = init_route, param = param_route )

  route_ice    <- Glacier_Disch(model = 1, inputData = cbind(ice_swe, ice_disch),
                                initCond = init_route_ice, param = param_route_ice  )

  # transfer function
  tf_soil <- round(
    UH(model = 1, Qg = route_module[ , 1], param = param_tf), 4  )

  tf_ice  <- round(
    UH(model = 1, Qg = route_ice[ , 1], param = param_tf), 4  )

  tf_out  <- tf_soil + tf_ice

  return( cbind(total = tf_out, soil = tf_soil, glacier = tf_ice) )


}# end fun