Appendix B: Prove for weights objection function

Contents

Appendix B: Prove for weights objection function#

The parameters and objective value for the calibration period for the different weight ratios are shown below.

General code#

#Loading packages
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

import numpy as np
from pathlib import Path
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr

import ewatercycle
import ewatercycle.models
import ewatercycle.forcing

from scipy.stats import qmc

#Loading shapefile of the catchment area of Mohembo
shapefile_path =  Path.home() / "BEP-beau/book/thesis_projects/BSc/2026_Q4_BeauBuijtenhuijs_CEG/Data/Shapefile" / "CatchmentArea_4326.shp"

#Defining start and end date for ERA5 data
start_date = "1970-01-01T00:00:00Z"
end_date = "2020-12-31T00:00:00Z"

#Creating path for ERA5 data
forcing_path_ERA5 = Path.home() / "BEP-beau/book/thesis_projects/BSc/2026_Q4_BeauBuijtenhuijs_CEG/Data/ERA5_HBV/own_shapefile_3"
forcing_path_ERA5.mkdir(exist_ok=True, parents=True)

#Downloading ERA5 data
#ERA5_forcing = ewatercycle.forcing.sources["LumpedMakkinkForcing"].generate(
#   dataset="ERA5",
#   start_time=start_date,
#   end_time=end_date,
#   shape=shapefile_path,
#   directory=forcing_path_ERA5,
#)

#Loading downloaded ERA5 data
load_location = forcing_path_ERA5 / "work" / "diagnostic" / "script"  
ERA5_forcing = ewatercycle.forcing.sources["LumpedMakkinkForcing"].load(directory=load_location)

#Defining start and end date for the calibration period
calibration_start_date = "1975-01-01T00:00:00Z"
calibration_end_date = "2010-12-31T00:00:00Z"

#Loading observed discharge data
data = pd.read_csv(Path.home() / "BEP-beau/book/thesis_projects/BSc/2026_Q4_BeauBuijtenhuijs_CEG/Data/mohembo_daily_water_discharge_data.csv", 
                   index_col='date', parse_dates=True, dayfirst=True)
data_daily = data.resample('D').interpolate()
data_daily.columns = ['Discharge (m$^3$/s)']
data_daily = data_daily[~data_daily.index.year.isin([1974, 2021])]

#The model returns data in mm/day, while observed data is in m^3/s
Area_km2 = 173696.852
def mmday_to_m3s(mmday_data, area):
    return (mmday_data * area) / 86.4

#Objective function
def objective(modelled, observed, start, end, w_kge, w_lognse):
    start = pd.to_datetime(start).tz_localize(None)
    end = pd.to_datetime(end).tz_localize(None)
    
    modelled.index = pd.to_datetime(modelled.index)
    observed.index = pd.to_datetime(observed.index)

    df = pd.concat([modelled.reindex(observed.index, method='ffill'), observed], axis=1, keys=['Modelled', 'Observed'])
    df = df.dropna()
    df = df[(df.index > start) & (df.index < end)]

    #KGE
    r = np.corrcoef(df['Observed'], df['Modelled'])[0, 1]
    beta = np.mean(df['Modelled']) / np.mean(df['Observed'])
    CV_modelled = np.std(df['Modelled']) / np.mean(df['Modelled'])
    CV_observed = np.std(df['Observed']) / np.mean(df['Observed'])                                 
    gamma = CV_modelled / CV_observed
    
    kge = 1 - np.sqrt((r - 1)**2 + (beta - 1)**2 + (gamma - 1)**2)

    #logNSE
    a = (np.log(df['Modelled']) - np.log(df['Observed']))**2
    a1 = a.sum() 
    b = (np.log(df['Observed']) - np.log(np.mean(df['Observed'])))**2
    b1 = b.sum()
    lognse = 1 - (a1/b1)

    score = w_kge * kge + w_lognse * lognse
    return score

#Making sets of parameters
N = 10 #2000 were used to calibrate the parameters, but this takes a long time
s_0 = np.array([0,  100,  0,  5,  0])

param_names = ["Imax", "Ce", "Sumax", "Beta", "Pmax", "Tlag", "Kf", "Ks", "Fm"]
param_mins = np.array([0, 0.2, 1000, 0.5, 0.001, 1, 0.0005, 0.00005, 0.00001])
param_maxs = np.array([8, 1.5, 3000, 4, 0.3, 90, 0.1, 0.002, 1])

sampler = qmc.LatinHypercube(d=len(param_names))
sample = sampler.random(n=N)
parameters = qmc.scale(sample, param_mins, param_maxs)

#Setting up the HBV model
ensemble = []
for counter in range(N): 
    ensemble.append(ewatercycle.models.HBVLocal(forcing=ERA5_forcing))
    config_file, _ = ensemble[counter].setup(parameters = parameters[counter],  initial_storage=s_0)
    ensemble[counter].initialize(config_file)

#Running the model with every generated parameters set
values = []
for ensembleMember in ensemble:
    Q_m = []
    time = []
    while ensembleMember.time < ensembleMember.end_time:
        ensembleMember.update()
        discharge_this_timestep = ensembleMember.get_value("Q")
        Q_m.append(discharge_this_timestep[0])
        time.append(ensembleMember.time_as_datetime)

    Q_m = mmday_to_m3s(np.array(Q_m), Area_km2)
    discharge_dataframe = pd.DataFrame({'model output': Q_m}, index=pd.to_datetime(time))

    #Calculating fit using the objective function
    #By changing the weights, results of different ratio can be computed
    obj_value = objective(discharge_dataframe['model output'], data_daily['Discharge (m$^3$/s)'], 
                          calibration_start_date, calibration_end_date, 0.3, 0.7)
    values.append(obj_value)
    
    del Q_m, time, discharge_dataframe, obj_value
    
for ensembleMember in ensemble:
    ensembleMember.finalize()

#Retrieving the best parameter set and best objective value for calibration period
best_index = np.argmax(values)
best_parameters = parameters[best_index]
best_value = values[best_index]

Weight KGE = 0.3, weight log-NSE = 0.7:#

Results
Imax	7.170 mm
Ce	0.589 (-)
Sumax	2653.448 mm
Beta	3.706 (-)
Pmax	0.214 mm/day
Tlag	85.687 days
Kf	0.078 day\(^{-1}\)
Ks	0.001 day\(^{-1}\)
Fm	0.468 (-)
Objective value	0.811

Weight KGE = 0.4, weight log-NSE = 0.6:#

Results
Imax	5.021 mm
Ce	0.594 (-)
Sumax	2012.102 mm
Beta	3.954 (-)
Pmax	0.219 mm/day
Tlag	87.626 days
Kf	0.067 day\(^{-1}\)
Ks	0.002 day\(^{-1}\)
Fm	0.731 (-)
Objective value	0.816

Weight KGE = 0.5, weight log-NSE = 0.5:#

Results
Imax	7.304 mm
Ce	0.546 (-)
Sumax	2683.134 mm
Beta	3.592 (-)
Pmax	0.182 mm/day
Tlag	82.451 days
Kf	0.042 day\(^{-1}\)
Ks	0.002 day\(^{-1}\)
Fm	0.845 (-)
Objective value	0.836

Weight KGE = 0.6, weight log-NSE = 0.4:#

Results
Imax	3.160 mm
Ce	0.518 (-)
Sumax	2141.544 mm
Beta	3.651 (-)
Pmax	0.256 mm/day
Tlag	67.745 days
Kf	0.082 day\(^{-1}\)
Ks	0.001 day\(^{-1}\)
Fm	0.061 (-)
Objective value	0.815

Weight KGE = 0.7, weight log-NSE = 0.3:#

Results
Imax	7.003 mm
Ce	0.412 (-)
Sumax	2248.938 mm
Beta	2.738 (-)
Pmax	0.185 mm/day
Tlag	22.48 days
Kf	0.015 day\(^{-1}\)
Ks	0.001 day\(^{-1}\)
Fm	0.372 (-)
Objective value	0.854

Weight KGE = 0.8, weight log-NSE = 0.2:#

Results
Imax	3.428 mm
Ce	0.211 (-)
Sumax	1180.128 mm
Beta	2.129 (-)
Pmax	0.276 mm/day
Tlag	15.057 days
Kf	0.014 day\(^{-1}\)
Ks	0.001 day\(^{-1}\)
Fm	0.669 (-)
Objective value	0.821