In the previous chapters we deal with ideal or synthetic cases. This purely academic examples, were though as a gentle (and pedagogical) introduction to the HBV.IANIGLA package. In this chapter we are going to face our first real world case: the simulation of a glacier surface mass balance.
The simulation of glacier mass balances is relevant in the Andes Cordillera, where these ice bodies have an important contribution to catchment discharge (Pellicciotti et al. 2008; Ragettli and Pellicciotti 2012; La Frenierre and Mark 2014; A. Ayala et al. 2016; Masiokas et al. 2016; Á. Ayala et al. 2020).
In mountain areas with scarce meteorological information, temperature index models are widely used to simulate snow and ice melting (Hock 2003; Konz and Seibert 2010; Finger et al. 2015; A. Ayala et al. 2017; Bravo et al. 2017). Since air temperature is the most readily available meteorological data in remote areas, the temperature index approach has been widely used in glaciological and hydrological modeling (Ohmura 2001). This package has been built with the SnowGlacier_HBV function, a module that uses this empirical approach to simulate snow, clean ice and debris covered melting.
In this section, we simulate the glacier mass balance for the Alerce glacier. Located on Monte Tronador (41.15º S ; 71.88º W), nearby the border between Argentina and Chile in the Andes of Northern Patagonia, Alerce is a medium size mountain glacier with an area of about 2.33 km2 that ranges between 1629 and 2358 m elevation showing a SE aspect (Ruiz et al. 2017; IANIGLA-ING 2018).
Since 2013 it is part of the monitoring network of the Inventario
Nacional de Glaciares (IANIGLA-ING
2010). Measurements are carried out following the glaciological
method for seasonal mass balance computation (Kaser, Fountain, and Jansson 2003).
Precipitation series from Puerto Montt (Dirección General de Aguas,
Chile) and air temperature series from Bariloche (Servicio
Meteorológico Nacional - Argentina) were used as meteorological
forcing to simulate the annual glacier mass balance
(data(alerce_data)).
When calibrating the model parameters, we considered as acceptable simulations those showing an annual mass balance in the range of MB ± 400 mm, where MB is the glacier annual mass balance.
library(HBV.IANIGLA)
## loading the data-set
data(alerce_data)
# now extract
meteo_data <- alerce_data$meteo_data # meteorological forcing series
mass_balance <- alerce_data$mass_balance # annual glacier mass balances
mb_dates <- alerce_data$mb_dates # fix seasonal dates
gl_topo <- alerce_data$topography # elevation bands
z_tair <- alerce_data$station_height[1] # topographic elevation air temp.
z_precip <- alerce_data$station_height[2] # topographic elevation of precipitation gauge
Despite the fact that we have a global surface mass balance estimation (derived from measurements), in order to take into consideration the topographic effect on it, we discretised our ice body in elevation bands. This mean that we have to build a semi-distributed glacier surface mass balance model.
# in this example we incorporate the precipitation and
# air temperature extrapolation functions into our model
## agruments description
# topograpphy: data frame with the elevation bands (gl_topo in this script).
# meteo: data frame with dates, air temperature and precipitation series.
# z_topo: numeric vector with air temperature and precipitation gauge elevation.
# param_tair: numeric vector with air temperature module parameters.
# param_precip: numeric vector with precipitation module parameters.
# param_ice: numeric vector with the glacier module parameters.
# init_ice: numeric value with initial snow water equivalent. Default value being 00 mm.
## output
# data frame with two columns: the date and the daily mass balance.
glacier_hbv <- function(topography,
meteo,
z_topo,
param_tair,
param_precip,
param_ice,
init_ice = 0){
n_it <- nrow(topography) # to get the number of elevation bands
n_dates <- nrow(meteo) # to get number of rows
precip_model <- tair_model <- matrix(NA_real_, nrow = n_dates, ncol = n_it)
glacier_model <- list()
for(i in 1:n_it){
# we first distribute the meteo forcing among the elevation bands
precip_model[ , i] <- Precip_model(model = 1,
inputData = meteo[ , "P(mm/d)"],
zmeteo = z_topo[2],
ztopo = topography[i , "mean"],
param = param_precip)
tair_model[ , i] <- Temp_model(model = 1,
inputData = meteo[ , "Tair(ºC)"],
zmeteo = z_topo[1],
ztopo = topography[i , "mean"],
param = param_tair)
# now we use the extrapolated values in the glacier module
glacier_model[[ i ]] <-
SnowGlacier_HBV(model = 1,
inputData = cbind( tair_model[ , i], precip_model[ , i]),
initCond = c(init_ice, 1, topography[i, "area_rel"]),
param = param_ice )
}
# we aggregate the mass balance series for the whole glacier
cum_mb <- lapply(X = 1:n_it, FUN = function(x){
out <- glacier_model[[x]][ , 7] * topography[x, "area_rel"]
})
cum_mb <- Reduce(f = `+`, x = cum_mb)
# return the column Cum = Psnow - Mtot (aggregated at glacier scale)
cum_out <- data.frame(date = meteo[ , "Date"],
"cum_mb(mm)" = cum_mb,
check.names = FALSE)
return(cum_out)
}When the measures and simulations are at different time-scale (annual vs daily) is always useful to construct an aggregation function (apart from the model).
# this function will use the glacier_hbv output data frame
# in order to get the annual mass balance
## arguments
# x: data frame resulting from glacier_hbv model
# start_date: string vector with starting dates.
# end_date: string vector with ending dates.
## output
# data frame with the annual mass balances
agg_mb <- function(x, start_date, end_date){
n_it <- length(start_date)
annual_mb <- c()
for(i in 1:n_it){
flag_start <- which(x[ , 1] == start_date[i])
flag_end <- which(x[ , 1] == end_date[i])
annual_mb[i] <- sum( x[flag_start:flag_end, 2] )
}
# build output data frame
out <- data.frame(start = start_date,
end = end_date,
"mb(mm we)" = annual_mb,
check.names = FALSE)
return(out)
}
# in this chunk of code I will also build the GOF function (you can also take a look
# at the hydroGOF package)
## arguments
# obs_upp: numeric vector with the acceptable upper limit
# obs_lwr: numeric vector with the acceptable lower limit
# sim: numeric vector with simulated mass balance
## output
# numeric value with the number of times that the simulation fits in between
# the lower and upper limits.
my_gof <- function(obs_upp, obs_lwr, sim){
n_it <- length(sim)
out <- 0
for(i in 1:n_it){
if(sim[i] >= obs_lwr[i] & sim[i] <= obs_upp[i]){
out <- 1 + out
} else {
out <- 0 + out
}
}
return(out)
}As in the previous chapters (3 and 4), we are going to sample the parameter space in order to get the acceptable parameter sets.
# air temperature model
tair_range <- rbind(
t_grad = c(-9.8, -2)
)
# precip model
precip_range <- rbind(
p_grad = c(5, 25)
)
# glacier module
glacier_range <- rbind(
sfcf = c(1, 2),
tr = c(0, 3),
tt = c(0, 3),
fm = c(1, 4),
fi = c(4, 8)
)
## we aggregate them in a matrix
param_range <-
rbind(
tair_range,
precip_range,
glacier_range
)In the next step we generate the random parameter sets,
# set the number of model runs that you want to try
n_run <- 10000
# build the matrix
n_it <- nrow(param_range)
param_sets <- matrix(NA_real_, nrow = n_run, ncol = n_it)
colnames(param_sets) <- rownames(param_range)
set.seed(123) # just for reproducibility
for(i in 1:n_it){
param_sets[ , i] <- runif(n = n_run,
min = param_range[i, 1],
max = param_range[i, 2]
)
}
head(param_sets)
#> t_grad p_grad sfcf tr tt fm fi
#> [1,] -7.556895 11.211834 1.991123 0.2636769 0.6044820 1.639739 5.692968
#> [2,] -3.651220 11.490401 1.302231 2.3388344 2.9345710 2.331002 7.141930
#> [3,] -6.609980 22.405084 1.433759 0.8435020 1.2381374 1.384493 6.407414
#> [4,] -2.912464 11.573476 1.160521 1.6490971 1.4259374 2.434075 4.327736
#> [5,] -2.464355 7.514025 1.823027 0.5838206 0.0531312 3.157230 5.180723
#> [6,] -9.444659 12.124429 1.208091 1.0183195 1.1042341 3.717423 5.354845Now we combine our functions and extract our best simulations,
# goodness of fit vector
gof <- c()
# make a loop
for(i in 1:n_run){
# run the model
glacier_sim <- glacier_hbv(topography = gl_topo,
meteo = meteo_data,
z_topo = c(z_tair, z_precip),
param_tair = param_sets[i, rownames(tair_range)],
param_precip = param_sets[i, rownames(precip_range) ],
param_ice = param_sets[i, rownames(glacier_range)] )
# aggregate the simulation
annual_mb <- agg_mb(x = glacier_sim,
start_date = as.Date( mb_dates$winter[-4] ),
end_date = as.Date( mb_dates$winter[-1] ) - 1 )
# compare the simulations with measurements
gof[i] <- my_gof(obs_upp = mass_balance$upp,
obs_lwr = mass_balance$lwr,
sim = annual_mb[ , 3])
rm(glacier_sim, annual_mb)
}
param_sets <- cbind(param_sets, gof)
# we apply a filter
param_subset <- subset(x = param_sets, subset = gof == 3)Once we have the subsetted our parameter matrix, we run the simulations to obtain a mean value (one per year).
# now we run the model again to get our simulations
n_it <- nrow(param_subset)
mb_sim <- matrix(NA_real_, nrow = 3, ncol = n_it)
for(i in 1:n_it){
glacier_sim <- glacier_hbv(topography = gl_topo,
meteo = meteo_data,
z_topo = c(z_tair, z_precip),
param_tair = param_subset[i, rownames(tair_range)],
param_precip = param_subset[i, rownames(precip_range) ],
param_ice = param_subset[i, rownames(glacier_range)] )
annual_mb <- agg_mb(x = glacier_sim,
start_date = as.Date( mb_dates$winter[-4] ),
end_date = as.Date( mb_dates$winter[-1] ) - 1 )
mb_sim[ , i] <- annual_mb[ , 3]
rm(i, glacier_sim, annual_mb)
}
# now we are going to make a data frame with the mean surface mass balance simulation
mean_sim <- cbind( mass_balance,
"mb_sim" = rowMeans(mb_sim) )
# make the plot
library(ggplot2)
g1 <-
ggplot(data = mean_sim, aes(x = year)) +
geom_pointrange(aes(y = `mb(mm we)`, ymin = `lwr`, color = 'obs',
ymax = `upp` ), size = 1, fill = "white", shape = 21) +
geom_point(aes(y = `mb_sim`, fill = 'sim'), shape = 23,
size = 3) +
geom_hline(yintercept = 0) +
scale_y_continuous(limits = c(-1500, 500), breaks = seq(-1500, 500, 250) ) +
scale_color_manual(name = '', values = c('obs' = 'blue') ) +
scale_fill_manual(name = '', values = c('sim' = 'red') ) +
ggtitle('') +
xlab('') + ylab('mb (mm we)') +
theme_minimal() +
theme(
title = element_text(color = "black", size = 12, face = "bold"),
axis.title.x = element_text(color = "black", size = 12, face = "bold"),
axis.title.y = element_text(color = "black", size = 12, face = "bold"),
legend.text = element_text(size = 11),
axis.text = element_text(size = 11))
g1Is your turn