5  Exercise 2: Implementation of the model

5.1 Assignment

The lumped, conceptual rainfall-runoff model implemented here is a simplified and modified version of the FLEX model as described in Fenicia et al. (2006) and Fenicia et al. (2008). The simplifications include not considering an interception reservoir and omitting lag functions for the routing. An overview of the model structure is given ADD FIGURE. The model considers three state variables, corresponding to the volume of water in three reservoirs:

1. The unsaturated soil reservoir \(S(t)\) [mm]
2. The fast reacting reservoir \(S_1(t)\) [mm]
3. The slow reacting reservoir \(S_2(t)\) [mm] , conceptually representing the groundwater system

The model is forced with precipitation \(P(t)\) and potential evaporation \(E_p(t)\) (see Section 3.3 for details on the forcing data). Note that mm is used as a unit here, as it allows to implement the model independently of the catchment area¹

Table 5.1: Model parameters together with their units, minimum and maximum value and a set of uncalibrated values. Note that \(S_{2,\mathrm{max}}\) can be used as the maximum capacity of both the slow and fast reacting reservoir.

Parameter	Units	Minimum	Maximum	Value
\(\lambda\)	-	0.5	1.5	0.65
\(S_{\mathrm{max}}\)	mm	87	430	350
\(b\)	-	0.5	1	0.86
\(\alpha\)	-	1	1.5	1.1
\(Q_{p,\mathrm{max}}\)	mm/d	0.75	75	34
\(\beta\)	-	0.05	0.1	0.06
\(\gamma\)	-	5	10	9.8
\(S_{2,\mathrm{max}}\)	mm	13	26	14
\(\kappa_2\)	mm/d	75	190	110
\(\kappa_1\)	1/d	0.26	0.52	0.32

By applying the conservation of mass to each reservoir, three ordinary differential equations can be derived to describe the change in storage in each reservoir over time.

\[ \frac{dS(t)}{dt} = P_{\mathrm{in}}(t) - E_{\mathrm{a}}(t) - Q_{\mathrm{p}}(t) \tag{5.1}\] \[ \frac{dS_1(t)}{dt} = P_1(t) - Q_1(t) + Q_{\mathrm{p}}(t) \tag{5.2}\] \[ \frac{dS_2(t)}{dt} = P_2(t) - Q_2(t) \tag{5.3}\]

In the equations above, a number of fluxes are calculated which depend on the state variables themselves. A first flux is the actual evaporation (\(E_{\mathrm{a}}\) [mm/d]), which can be calculated from the potential evaporation \(E_{\mathrm{p}}\) as follows:

\[ E_{\mathrm{a}}(t) = \frac{1}{\lambda}\frac{S(t)}{S_{\mathrm{max}}}E_{\mathrm{p}}(t) \tag{5.4}\]

where \(\lambda\) [-] is a dimensionless model parameter and \(S_{\mathrm{max}}\) [mm] is the storage capacity of the unsaturated soil reservoir \(S(t)\). The infiltration in the soil reservoir \(P_{\mathrm{in}}\) [mm/d] is calculated using:

\[ P_{\mathrm{in}}(t) = \left(1 - \frac{S(t)}{S_{\mathrm{max}}}\right)^b P(t) \tag{5.5}\]

where \(b\) is again a dimensionless model parameter. Subsequently, the excess precipitation \(P_{\mathrm{exc}}\) [mm/d] can be estimated:

\[ P_{\mathrm{exc}}(t) = P(t) - P_{\mathrm{in}}(t) \tag{5.6}\]

The amount of water percolating to the groundwater storage \(Q_{\mathrm{p}}\) [mm/d] is:

\[ Q_{\mathrm{p}}(t)=Q_{\mathrm{p,max}}\left(1-e^{-\beta\frac{S(t)}{S_{\mathrm{max}}}}\right) \tag{5.7}\]

with \(Q_{\mathrm{p,max}}\) [mm/d] the maximum percolation volume per unit of time and \(\beta\) [-] a dimensionless parameter.

\(P_{\mathrm{exc}}\) is again partitioned between \(S_2\) and \(S_1\): \[ P_{\mathrm{exc}}(t) = P_1(t) + P_2(t) \tag{5.8}\] The part of the effective precipitation reaching \(S_2\) is:

\[ P_2(t) = \alpha\frac{S(t)}{S_{\mathrm{max}}}P_{\mathrm{eff}}(t) \]

with \(\alpha\) [-] a dimensionless model parameter. \(P_1\) is then calculated as the remaining part of the effective precipitation using Equation 5.8.

\(S_2\) is considered as a non-linear reservoir, so that its outflow \(Q_2\) [mm/d] can be calculated using:

\[ Q_2(t) = \kappa_2\left(\frac{S_2(t)}{S_{2,\mathrm{max}}}\right)^\gamma \tag{5.9}\]

with \(S_{2,\mathrm{max}}\) [mm] the storage capacity of the fast reacting reservoir, \(\kappa_2\) [mm/d] the maximum outflow rate, and \(\gamma\) [-] a dimensionless model parameter.

\(S_1\) on the other hand is considered as a linear reservoir, so the outflow \(Q_1\) [mm/d] is given by: \[ Q_1(t) = \kappa_1 S_1(t) \tag{5.10}\]

where \(\kappa_1\) [1/d] the reciprocal of the residence time. The total discharge \(Q\) [mm/d] is then calculated as the sum of the outflow from the fast and slow reacting reservoir: \[ Q(t) = Q_1(t) + Q_2(t) \tag{5.11}\]

Equation 5.1, Equation 5.2, and Equation 5.3 can be denoted as a general non-linear, continuous-time state-space model: \[ \frac{d\mathbf{x}}{dt} = \mathcal{f}(\mathbf{x}, \mathbf{{f}}, \mathbf{p}) \]

with \(f\) a non-linear function, \(\mathbf{x} = [S, S_1, S_2]^T\) the state vector, \(\mathbf{f} = [P, E_p]^T\) the forcing vector, and \(\mathbf{p} = [\lambda, S_{\mathrm{max}}, b, \alpha, Q_{\mathrm{p,max}}, \beta, \gamma, S_{2,\mathrm{max}}, \kappa_2, \kappa_1]^T\) the parameter vector. For this exercise, the model will be implemented in discrete time using a simple forward Euler scheme²: \[ \mathbf{x}(t+\Delta t) = \mathbf{x}(t) + \mathcal{f}(\mathbf{x}(t), \mathbf{{f}}(t), \mathbf{p})\Delta t \] The time step \(\Delta t\) is set to 1 day, which is the same as the time step of the forcing data.

With all of the information above, implement the model based on Equation 5.1 to Equation 5.11. As a start, run it with the uncalibrated parameter values given in Table 5.1. Compare the simulated discharge with the observed discharge (see Section 3.3). Discuss the performance using relevant figures and metrics.

Tip

Some general tips and recommendations:

• Implement the model in a separate file scripts/model.py as a function. You need to be able to reuse the model in the following exercises.
• Make sure you can run the model one time step at a time. A for-loop can then be used to run the model for the entire time period.
• Take into account the physical limitations of the state variables.
• Take into account the impact of initialisation of the state variables

5.2 Results and discussion

Fenicia, Fabrizio, Hubert H. G. Savenije, Patrick Matgen, and Laurent Pfister. 2008. “Understanding Catchment Behavior Through Stepwise Model Concept Improvement.” Water Resources Research 44 (1). https://doi.org/10.1029/2006WR005563.

Fenicia, F., H. H. G. Savenije, P. Matgen, and L. Pfister. 2006. “Is the Groundwater Reservoir Linear? Learning from Data in Hydrological Modelling.” Hydrology and Earth System Sciences 10 (1): 139–50. https://doi.org/10.5194/hess-10-139-2006.

---

1. \(1 \mathrm{ mm} = 1 \cdot 10^{-3} \mathrm{ m}^3 / \mathrm{m}^2\)

2. For more background info on numerical methods for ordinary differential equation, see your bachelor course on differential equations