Inferring X4 in GR4J without using the model
<!DOCTYPE html>
Inferring X4 in GR4J without using the model¶
Introduction¶
In this notebook I will show that it is not necessary to run the GR4J hydrological model in order to infer its X4 parameter, as it is usually done when performing brute-force calibration. Instead, the lag between rainfall events and the resulting streamflow peaks is statistically computed and used as a proxy for the value of X4. If calibration is used to infer the other parameters of the model, this approach allows to remove one dimension in the parameter space and consequently to save some CPU. Also, by fixing a parameter (hopefully) to its true value, it decreases the ability of the model to compensate errors with other parameters, thus reducing the chances of equifinality and increasing the robustness of the model.
Material¶
This notebook is available on GitHub. It runs on Python3 with a number of libraries that can be installed through the conda package manager. I will use time series of precipitation (P), potential evapotranspiration (E) and discharge (Q), which cover the 1948-1988 time range at a daily time step.
import random
from tqdm import tqdm
import numpy as np
from numba import jit
from scipy.interpolate import interp1d
from scipy.optimize import fmin
import pandas as pd
from pandas import DataFrame
from matplotlib import pyplot as plt
from IPython.display import Image
%matplotlib inline
def plot_peq(peq):
fig, ax1 = plt.subplots(figsize=(15, 5))
ax1.set_ylim([peq.p.max() * 2, 0])
ax1.plot(peq.p, color='b', alpha=0.5, label='Rainfall')
ax1.set_ylabel('Rainfall (mm/day)', color='b')
for tl in ax1.get_yticklabels():
tl.set_color('b')
ax2 = ax1.twinx()
ax2.set_ylim([0, peq.q.max() * 2])
ax2.plot(peq.e, color='r', alpha=0.3, label='PET')
ax2.plot(peq.q, color='g', label='Discharge')
ax2.set_ylabel('Discharge (mm/day)', color='g')
for tl in ax2.get_yticklabels():
tl.set_color('g')
ax1.legend(loc='center left')
ax2.legend(loc='center right')
plt.show()
peq_full = pd.read_pickle('peq.pkl')
plot_peq(peq_full.loc['1985':'1987'])
The code for the GR4J model (represented in the figure below) is provided here in pure Python, and Just-In-Time compiled using Numba for performance reasons.
Image(filename='gr4j.png')
@jit
def run_gr4j(x, p, e, q, s, uh1_array, uh2_array, l, m):
for t in range(p.size):
if p[t] > e[t]:
pn = p[t] - e[t]
en = 0.
tmp = s[0] / x[0]
ps = x[0] * (1. - tmp * tmp) * np.tanh(pn / x[0]) / (1. + tmp * np.tanh(pn / x[0]))
s[0] += ps
elif p[t] < e[t]:
ps = 0.
pn = 0.
en = e[t] - p[t]
tmp = s[0] / x[0]
es = s[0] * (2. - tmp) * np.tanh(en / x[0]) / (1. + (1. - tmp) * np.tanh(en / x[0]))
tmp = s[0] - es
if tmp > 0.:
s[0] = tmp
else:
s[0] = 0.
else:
pn = 0.
en = 0.
ps = 0.
tmp = (4. * s[0] / (9. * x[0]))
perc = s[0] * (1. - (1. + tmp * tmp * tmp * tmp) ** (-1. / 4.))
s[0] -= perc
pr_0 = perc + pn - ps
q9 = 0.
q1 = 0.
for i in range(m):
if i == 0:
pr_i = pr_0
else:
pr_i = s[2 + i - 1]
if i < l:
q9 += uh1_array[i] * pr_i;
q1 += uh2_array[i] * pr_i;
q9 *= 0.9
q1 *= 0.1
f = x[1] * ((s[1] / x[2]) ** (7. / 2.))
tmp = s[1] + q9 + f
if tmp > 0.:
s[1] = tmp
else:
s[1] = 0.
tmp = s[1] / x[2]
qr = s[1] * (1. - ((1. + tmp * tmp * tmp * tmp) ** (-1. / 4.)))
s[1] -= qr
tmp = q1 + f
if tmp > 0.:
qd = tmp
else:
qd = 0.
q[t] = qr + qd
if s.size > 2:
s[3:] = s[2:-1]
s[2] = pr_0
class gr4j:
def sh1(self, t):
if t == 0:
res = 0.
elif t < self.x[3]:
res = (float(t) / self.x[3]) ** (5. / 2.)
else:
res = 1.
return res
def sh2(self, t):
if t == 0:
res = 0.
elif t < self.x[3]:
res = 0.5 * ((float(t) / self.x[3]) ** (5. / 2.))
elif t < 2. * self.x[3]:
res = 1. - 0.5 * ((2. - float(t) / self.x[3]) ** (5. / 2.))
else:
res = 1.
return res
def uh1(self, j):
return self.sh1(j) - self.sh1(j - 1)
def uh2(self, j):
return self.sh2(j) - self.sh2(j - 1)
def __init__(self, x):
self.x = np.array(x)
self.s = np.empty(2 + int(2. * self.x[3]))
self.s[0] = self.x[0] / 2.
self.s[1] = self.x[2] / 2.
self.s[2:] = 0.
self.l = int(self.x[3]) + 1
self.m = int(2. * self.x[3]) + 1
self.uh1_array = np.empty(self.l)
self.uh2_array = np.empty(self.m)
for i in range(self.m):
if i < self.l:
self.uh1_array[i] = self.uh1(i + 1)
self.uh2_array[i] = self.uh2(i + 1)
def run(self, pe):
q = np.empty_like(pe[0])
run_gr4j(self.x, pe[0], pe[1], q, self.s, self.uh1_array, self.uh2_array, self.l, self.m)
return [q]
Method¶
We will start by fitting the model to the observed data with a brute-force calibration.
def calibration(x, in_obs, out_obs, warmup_period, crit_func, model, x_range, x_fix):
this_x = []
for i in range(len(x_range)):
if x_fix[i] is None:
if x[i] < x_range[i][0]:
return np.inf
if x[i] > x_range[i][1]:
return np.inf
this_x.append(x[i])
else:
this_x.append(x_fix[i])
data_nb = out_obs[0].size
q_mod = model(this_x)
out_sim = q_mod.run(in_obs)
error = crit_func(out_obs[0][warmup_period:], out_sim[0][warmup_period:])
return error
def nse(x_obs, x_est):
return 1. - (np.sum(np.square(x_obs - x_est)) / np.sum(np.square(x_obs - np.mean(x_obs))))
def nse_min(x_obs, x_est):
return 1. - nse(x_obs, x_est)
peq = peq_full
x_range = [[0., np.inf], [-np.inf, np.inf], [0., np.inf], [0., np.inf]]
x_fix = [None, None, None, None]
x0 = [100, 0, 100, 1]
warmup_period = 365
x = fmin(calibration, x0, args=([peq.p.values, peq.e.values], [peq.q.values], warmup_period, nse_min, gr4j, x_range, x_fix), xtol=0.01, ftol=0.01, disp=False)
x4_real = x[3]
q_mod = gr4j(x)
peq['q_sim'] = q_mod.run([peq.p.values, peq.e.values])[0]
ax = peq.loc['1985':'1987', ['q', 'q_sim']].plot(figsize=(15, 5))
ax.set_ylabel('mm/day')
print('X = ' + str(x))
print('NSE = ' + str(nse(peq.q.values, peq.q_sim)))
As can be seen from the NSE score (0.87) and the hydrographs in the figure above, the simulated streamflow fits the observed streamflow quite well. So from this model calibration, the inferred value is X4 = 3.5.
Because X4 is the parameter that takes into account the lag time between precipitation and streamflow, it must be possible to measure it using only these two series. In fact, the cross-correlation does just that. We will limit the computation from a lag time of 0 (it is the precipitation which comes first, and then the streamflow, so a negative lag time has no meaning), to a lag time of several days (because we know the propagation of the precipitation by runoff into the streamflow at the outlet of the basin is a matter of days at most, depending on the size of the basin of course).
@jit
def correlate(p, q, width):
corr = np.zeros(width)
for i in range(width):
size = len(p) - i
for j in range(size):
corr[i] += p[j] * q[i + j]
corr[i] /= size
return corr
peq = peq_full
day_nb = 10
corr = correlate(peq.p.values, peq.q.values, day_nb)
df1 = DataFrame()
df1['Lag time (days)'] = np.arange(day_nb)
df1['Cross-correlation'] = corr
df1 = df1.set_index('Lag time (days)')
df1.plot(figsize=(15, 5))
Note that the time step of the cross-correlation is the same as the one of the series which were used to compute it (here, one day), and we know that X4 can take real values (i.e. non-integer values). This problem appart, it seems that the lag time corresponding to the maximum value of the cross-correlation is linked to the value of X4. Without demonstrating it mathematically, we will show that adding 1 to this value leads to a good approximation of X4.
In order to infer not only integer values of X4, which would be too gross of an approximation, we will use spline functions to interpolate the cross-correlation.
f = interp1d(df1.index, df1['Cross-correlation'], kind='cubic')
df2 = DataFrame()
precision = 0.1
df2['Lag time (days)'] = np.arange(0, day_nb - 1, precision)
df2['Interpolated cross-correlation'] = f(df2['Lag time (days)'])
df2 = df2.set_index('Lag time (days)')
ax = df1.plot(figsize=(15, 5))
df2.plot(ax=ax, figsize=(15, 5))
xmax = df2['Interpolated cross-correlation'].argmax()
ymax = df2['Interpolated cross-correlation'].loc[xmax]
ax.scatter(xmax, ymax, color='r', s=50)
x4_inferred = xmax + 1
ax.annotate('X4 inferred = ' + str(xmax) + ' + 1 = ' + str(x4_inferred) + '\n(should be ' + str(round(x4_real, 1)) + ')', xy=(xmax, ymax), xytext=(10, 10), textcoords='offset points')
plt.xticks(list(plt.xticks()[0]) + [xmax])
ylim = ax.get_ylim()
ax.plot([xmax, xmax], ylim)
plt.ylim(ylim)
plt.show()
The inferred value of X4 using cross-correlation (3.8) is not so different than the one using brute-force calibration (3.5).
In order to check if this approach works for various values of X4, we will use the GR4J model to generate a broad range of streamflows. We will use the same precipitation and potential evapotranspiration time series, pick some values at random for X4, but also for X1, X2 and X3, simulate a corresponding streamflow time series, and infer X4 using the cross-correlation between the precipitation and the simulated streamflow. This way, we can measure the difference between the inferred X4 value and the "true" X4 value.
peq = peq_full
day_nb = 50
results = {'x1': [], 'x2': [], 'x3': [], 'x4': [], 'x4_inferred': [], 'nse1': []}
for i in tqdm(range(1000)):
x1 = random.uniform(100, 500)
x2 = random.uniform(-1, 1)
x3 = random.uniform(10, 100)
x4 = random.uniform(0.1, day_nb / 2)
q_mod = gr4j([x1, x2, x3, x4])
q_sim = q_mod.run([peq.p.values, peq.e.values])[0]
corr = correlate(peq.p.values, q_sim, day_nb)
f = interp1d(np.arange(day_nb), corr, kind='cubic')
corr = f(np.arange(0, day_nb - 1, precision))
x4_inferred = corr.argmax() * precision + 1
q_mod_inferred = gr4j([x1, x2, x3, x4_inferred])
q_sim_inferred = q_mod_inferred.run([peq.p.values, peq.e.values])[0]
results['x1'].append(x1)
results['x2'].append(x2)
results['x3'].append(x3)
results['x4'].append(x4)
results['x4_inferred'].append(x4_inferred)
results['nse1'].append(nse(q_sim[warmup_period:], q_sim_inferred[warmup_period:]))
results= DataFrame(results)
plt.close('all')
plt.figure(figsize=(10, 10))
plt.plot([0.1, day_nb / 2], [0.1, day_nb / 2], label='y = x')
plt.scatter(results['x4'], results['x4_inferred'], s=5, color='r', label='X4 inferred vs simulated')
plt.xlim(0, day_nb / 2)
plt.ylim(0, day_nb / 2)
plt.xlabel('X4 simulated')
plt.ylabel('X4 inferred')
plt.legend(loc='upper left')
plt.show()
The figure above shows that there is a good match between the inferred value of X4 and the true value. There is, however, a problem for X4 < 1. The inferred value is always X4 = 1 in this case. This is because with a time step of one day, there is no way to say when the streamflow peak resulting from a precipitation event on a given day will appear within that day. We can just say that it will happen the same day. But does it even matter? To see if it does, we will compute the NSE of the simulated streamflow using the model where X4 was inferred, compared to the simulated streamflow using the reference model (the one with the true X4).
plt.close('all')
plt.figure(figsize=(15, 5))
nse_vs_x4 = results.loc[results['x4'] <= 25].set_index(['x4']).sort_index()
plt.plot([1, 1], [-2, 2], '--', color='gray')
plt.plot(nse_vs_x4.index, nse_vs_x4['nse1'], color='b')
plt.ylim(-0.1, 1.1)
plt.xlabel('X4')
plt.ylabel('NSE')
plt.xticks(list(plt.xticks()[0]) + [1])
plt.title('Performances of the model using the inferred X4 value as a function of the true X4 value')
plt.show()
As we can see, there is actually no problem: the model is not sensitive to a value of X4 < 1. It is the same as X4 = 1 (NSE = 1). There is however a performance degradation for X4 ~ 1 (and X4 >= 1). This is where the errors in the inferred value of X4 are relatively bigger, compared to the absolute value of X4. But this is still acceptable (NSE ~ 0.85). Note that we are comparing the inferred valude of X4 to its "true" value, which is a purely theoretical situation that will never happen in real life. With a measured (and not simulated) streamflow, the value of X4 obtained with brute-force calibration would not be "more true" than the one obtained with the cross-correlation.
Still, it is interesting to see if the model can compensate the errors in X4 with its other parameters. We didn't let the choice to GR4J with the values of X1, X2 and X3. Let us now calibrate them, fixing X4 to the value inferred by cross-correlation.
peq = peq_full
day_nb = 50
for i in tqdm(range(len(results))):
x1 = results.loc[i, 'x1']
x2 = results.loc[i, 'x2']
x3 = results.loc[i, 'x3']
x4 = results.loc[i, 'x4']
q_mod = gr4j([x1, x2, x3, x4])
q_sim = q_mod.run([peq.p.values, peq.e.values])[0]
x4_inferred = results.loc[i, 'x4_inferred']
x_range = [[1., np.inf], [-np.inf, np.inf], [0.1, np.inf], [0.1, np.inf]]
x_fix = [None, None, None, x4_inferred]
x0 = [100, 0, 100]
x_inferred = fmin(calibration, x0, args=([peq.p.values, peq.e.values], [q_sim], warmup_period, nse_min, gr4j, x_range, x_fix), xtol=0.01, ftol=0.01, disp=False)
q_mod_inferred = gr4j(list(x_inferred) + [x4_inferred])
q_sim_inferred = q_mod_inferred.run([peq.p.values, peq.e.values])[0]
results.loc[i, 'x1_inferred'] = x_inferred[0]
results.loc[i, 'x2_inferred'] = x_inferred[1]
results.loc[i, 'x3_inferred'] = x_inferred[2]
results.loc[i, 'nse2'] = nse(q_sim[warmup_period:], q_sim_inferred[warmup_period:])
plt.close('all')
plt.figure(figsize=(15, 5))
nse_vs_x4 = results.loc[results['x4'] <= 25].set_index(['x4']).sort_index()
plt.plot([1, 1], [-1, 2], '--', color='gray')
plt.plot(nse_vs_x4.index, nse_vs_x4['nse1'], color='b', label='without optimization of X1, X2, X3', alpha=0.5)
plt.plot(nse_vs_x4.index, nse_vs_x4['nse2'], color='r', label='with optimization of X1, X2, X3')
plt.plot(nse_vs_x4.index, nse_vs_x4['nse2'] - nse_vs_x4['nse1'], color='g', label='NSE gain with optimization of X1, X2, X3')
plt.ylim(-0.1, 1.1)
plt.xlabel('X4')
plt.ylabel('NSE')
plt.xticks(list(plt.xticks()[0]) + [1])
plt.title('Performances of the model using inferred X4 as a function of true X4')
plt.legend(loc='center right')
plt.show()
We can see that there is not much gain, which is reassuring in a way: it means that in the GR4J parameter space, X4 is orthogonal to the other parameters, and thus, useful and needed.
Conclusion¶
In this notebook, I showed that it is not necessary to infer all the parameters of the GR4J model using brute-force calibration. In particular, X4 can be inferred using the cross-correlation between the precipitation and the streamflow time series. This leads to an easier calibration of the other parameters, because one dimension can be removed from the parameter space.
According to the GR4J model conceptor, Claude Michel, who didn't want to give a physical meaning to its parameters (and this is the reason why he called them X), we should not try to interpret the concepts behind the model and make them match the reality. But this is exactly what we have done here: X4 is definitely linked to the speed at which the rainfall runs off to the outlet of the basin. And I cannot imagine it could be otherwise, even for the other parameters. To me, with such a small number of parameters, it is not possible that a conceptual model simulates so well the streamflows we measure on the ground. Its concepts must also mimic the physical processes.