import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

HBV

Let us say we want to do research into the future of the Speyside distileries using HBV. We see that, using the caravan map, Scotland is part of Caravan, and some of Speyside is in there. So we will use that region to do our research on, using the HBV model.

# General python
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

# Niceties
from rich import print

# General eWaterCycle
import ewatercycle
import ewatercycle.models
import ewatercycle.forcing

# Region selection using Caravan
"""Avon at Delnashaugh, in Scotland"""
caravan_id = "camelsgb_8004"

# Time periods for experiment, note they are hydrological years
experiment_start_date = "2010-08-01T00:00:00Z"
experiment_end_date = "2014-08-31T00:00:00Z"

Now we will get the forcing data for our selected region.

# Location of caravan forcing
forcing_path_caravan = Path.home() / "forcing" / caravan_id / "caravan"
forcing_path_caravan.mkdir(exist_ok=True, parents=True)

# Generate the caravan forcing
caravan_forcing_object = ewatercycle.forcing.sources['CaravanForcing'].generate(
    start_time=experiment_start_date,
    end_time=experiment_end_date,
    directory=forcing_path_caravan,
    basin_id=caravan_id,
)
print(caravan_forcing_object)

# Note that you want to do this in a separate notebook, then load it in when you want to use it!!!
caravan_forcing_object = ewatercycle.forcing.sources["CaravanForcing"].load(directory=forcing_path_caravan)

CaravanForcing(
    start_time='2010-08-01T00:00:00Z',
    end_time='2014-08-31T00:00:00Z',
    directory=PosixPath('/home/mmelotto/forcing/camelsgb_8004/caravan'),
    shape=PosixPath('/home/mmelotto/forcing/camelsgb_8004/caravan/camelsgb_8004.shp'),
    filenames={
        'pr': 'camelsgb_8004_2010-08-01_2014-08-31_pr.nc',
        'evspsblpot': 'camelsgb_8004_2010-08-01_2014-08-31_evspsblpot.nc',
        'tasmax': 'camelsgb_8004_2010-08-01_2014-08-31_tasmax.nc',
        'tas': 'camelsgb_8004_2010-08-01_2014-08-31_tas.nc',
        'tasmin': 'camelsgb_8004_2010-08-01_2014-08-31_tasmin.nc',
        'Q': 'camelsgb_8004_2010-08-01_2014-08-31_Q.nc'
    }
)

The model run

# Setting the HBV parameters, these are by no means optimized
hbv_parameters = [7.085,  # Imax
                  0.837,  # Ce
                  76.373,  # Sumax
                  1.112,  # Beta
                  0.245,  # Pmax
                  7.801,  # Tlag
                  0.096,  # Kf
                  0.003,  # Ks
                  0.226  # FM
                  ]

# Set initial state values
s_0 = np.array([0,  # Si
                100,  # Su
                0,  # Sf
                5,  # Ss
                0  # Sp
                ])

# Create model object, notice the forcing object.
HBV_model = ewatercycle.models.HBV(forcing=caravan_forcing_object)

# Create config file in model.setup()
config_file, _ = HBV_model.setup(parameters=hbv_parameters, initial_storage=s_0)

# Initialize model
HBV_model.initialize(config_file)

# Run model, capture calculated discharge and timestamps
Q_m = []
time = []

while HBV_model.time < HBV_model.end_time:
    HBV_model.update()
    Q_m.append(HBV_model.get_value("Q")[0])
    time.append(pd.Timestamp(HBV_model.time_as_datetime))

# Finalize model (shuts down container, frees memory)
HBV_model.finalize()

Look at the results

# Make a pandas series of the modelled discharge
model_output = pd.Series(data=Q_m, name="Modelled_discharge", index=time)

# Load the observations from the caravan object
caravan_discharge_observation = xr.open_mfdataset([caravan_forcing_object['Q']])

# Make a plot of both the observations and modelled discharge
model_output.plot()
caravan_discharge_observation["Q"].plot(label="Observed discharge")
plt.legend()
plt.ylabel("Discharge (mm/d)")

Text(0, 0.5, 'Discharge (mm/d)')

# We want to also be able to use the output of this model run in different analyses. Therefore, we save it as a NetCDF file.
xr_model_output = model_output.to_xarray()

xr_model_output.attrs['units'] = 'mm/d'

# Save the xarray Dataset to a NetCDF file
xr_model_output.to_netcdf(f'~/forcing/{caravan_id}/river_discharge_data.nc')

Your research ahead

Now you have this modelled discharge, you can perform your own research. Some examples, keeping in mind you have a calibrated model:

• One case could be that you look at future climate data from CMIP to assess the need for:

  1. More distilleries

  2. Less whisky production in the future

  3. What is the best way forward, regarding climate scenarios, for more whisky

• Assess the current water intake by the distilleries.

• Does snow-melt affect whisky production.

• Can HBV predict seasonal low-flow periods that might coincide with peak whisky tourism demand?