Python and C package to solve the reservoir differential equation using a piecewise quadratic interpolation following the QuaSoARe method.
This package implements the Quadratic Solution of the Approximate Reservoir Equation (QuaSoARe) method described in the following paper: Lerat, J. (2025), Technical note: Quadratic Solution of the Approximate Reservoir Equation (QuaSoARe), Hydrol. Earth Syst. Sci., 29, 2003–2021, https://doi.org/10.5194/hess-29-2003-2025, 2025.
- Create a suitable python environment. We recommend using miniconda combined with the environment specification provided in the env_pyquasoare.yml file in this repository.
- Git clone this repository and run
pip install .
Solution of the production store from the GR4J daily rainfall-runoff model using QuaSoAre:
from pathlib import Path
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pyquasoare import approx, models
from hydrodiy.io import csv
# Package root path (might need modification)
froot = Path(__file__).parent.parent
# The production store of the GR4J model is characterised by
# the following differential equation:
# dS / dt = P (1 - [S/X1]**2) - E S/X1 (2 - S/X1) - a (S/X1)**5
# where S is the store volume (mm), X1 is the store capacity,
# P and E are the rainfall and evapotranspiration (mm/day) and
# a is constant set to 2.25**4/4 (=6.407).
#
# If we introduce the following variables:
# p = P/X1
# e = E/X1
# u = S/X1
# the previous equation becomes:
# du / dt = p (1 - x**2) - e x (2 -x) - a x**5
# this equation has 3 fluxes:
# * rainfall stored in store = p (1 - x**2)
# * actual evapotranspiration = - e x (2-x)
# * percolation = -a x**5
X1 = 400
fluxes = [
lambda x: 1 - x**2,
lambda x: -x*(2-x),
lambda x: -2.25**4/4*x
]
# We are now solving this differential equation with QuaSoARe:
# Definition of interpolation points
nalphas = 20
alphas = np.linspace(0., 1.2, nalphas)
# Quadratic piecewise interpolation of the flux functions
amat, bmat, cmat, cst = approx.quad_coefficient_matrix(fluxes, alphas)
# Creating random rainfall and PET data
nval = 1000
rain = np.maximum(np.random.exponential(8, size=nval) - 2, 0)
evap = 2 + 2 * (np.sin(np.arange(nval)/365.25 * 2 * math.pi) + 1)/2
# GR4J applies an interception function. This
# leads to
rain_intercept = np.maximum(rain - evap, 0.)
evap_intercept = np.maximum(evap - rain, 0.)
# The scalings indicated below correspond to variables
# 'p' and 'e' of the previous equation:
scalings = np.column_stack([rain_intercept/X1,
evap_intercept/X1,
np.ones(nval)])
# Run the model using QuaSoare
s0 = 1./2
niter, s1, fx = models.quad_model(alphas, scalings, \
amat, bmat, cmat, s0, 1.)
# All fluxes computed by QuaSoARe needs to be rescaled
# X1 because the equation was solved for variables
# divided by X1 (see equations above)
sims = np.column_stack([s1*X1, fx[:, 0]*X1, \
-fx[:, 1]*X1, -fx[:, 2]*X1])
# Plot results
plt.close("all")
fig, axs = plt.subplots(nrows=4, figsize=(15, 10), layout="constrained")
for iax, ax in enumerate(axs):
ax.plot(time, sims[:, iax])
plt.show()All results presented in the QuaSoARe paper can be generated by running the python script models_run.py. This script applies QuaSoARe to a set of test cases defined by the script argument '-t' which varies from 0 to 29. The test cases includes application of QuaSoARe to
- 6 catchments locaed in Eastern Australia,
- 5 hydrological models
To run all cases, the script needs to be launched within a loop as follows (assuming Linux/Mac OS bash script):
for taskid in {0..29}; do
python scripts/quasoare_paper_2024/models_run.py -t taskid
done Once the results are generated, the figures of the paper can be generated using the script figures_generate_all.py.
This project is licensed under the MIT License, which allows for free use, modification, and distribution of the code under the terms of the license.
For proper citation of this project, please refer to the CITATION.cff file, which provides guidance on how to cite the software and relevant publications.