Report

Contents

Report#

Part 1#

1.1 (0.5 points)

Classify the model equation based on the categories listed under chapter 2.1 of the MUDE textbook (paste your answer from task 1.1).

The model equation is a first-order, non-linear ODE. Two terms make it non-linear: \(S^\alpha\) and \(\frac{S_{i+1}}{S_{i+1} + S_c}\)

Rubric:

  • 0.5 if all three aspects are correct

  • 0.25 points if 2/3 aspects are correct

Part 2#

2.1 (0.5 points)

Write down the discretized equation for \(S_{i+1}\) based on the implicit Euler scheme (paste your result from task 2.1). Explain why an iterative scheme is needed to solve for \(S_{i+1}\).

\[S_{i+1} = S_i + \Delta t \left(P(t_{i+1}) - \text{ET}_{\text{pot}}(t_{i+1}) \frac{S_{i+1}}{S_{i+1} + S_c} - k {S_{i+1}}^\alpha\right)\]

An iterative scheme is needed because the equation is non-linear, and \(S_{i+1}\) appears on both sides of the equation. In this case, we can not simply bring all instances of \(S_{i+1}\) to one side of the equation – a closed form solution is not generally available.

Rubric:

  • 0.25 points for the correct equation

  • additional 0.25 points for an explanation why an iterative scheme is needed

2.2 (1 point)

Write down the function \(g(z)\) whose root is solved with Newton–Raphson, and specify what \(z\) represents in this context (use your solution from task 2.2).

For this model, the value \(z\) that we search is \(S_{i+1}\), the state at the new time step.

We want to find \(S_{i+1}\) where the following expression is zero, thus, fulfilling the equation for the implicit Euler update.

\[g(z) = q(S_{i+1}) = S_{i+1} - S_i - \Delta t f(S_{i+1}, t_{i+1})\]

Plugging in the ODE for our model yields:

\[g(z) = q(S_{i+1}) = S_{i+1} - S_i - \Delta t \left(P(t_{i+1}) - \text{ET}_{\text{pot}}(t_{i+1}) \frac{S_{i+1}}{S_{i+1} + S_c} - k {S_{i+1}}^\alpha\right)\]

Rubric:

  • 0.5 points for the correct expression (only 0.25 if the equation was copied from the book without plugging in the ODE)

  • additional 0.5 points for relating \(z\) to the storage at the new time point \(S_{i+1}\)

2.3 (0.5 point)

Derive and present the expression for \(g'(z)\) used in the Newton–Raphson update for the ODE of the conceptual hydrological model (use your answer from task 2.3).

\[\begin{split} \begin{aligned} g^\prime(z) = q^\prime(S_{i+1}) &= 1 - \Delta t \left(0 -\text{ET}_{\text{pot}}(t_{i+1}) \frac{S_c}{(S_{i+1} + S_c)^2} - \alpha k {S_{i+1}}^{\alpha-1}\right)\\ &= 1 + \Delta t \left(\text{ET}_{\text{pot}}(t_{i+1}) \frac{S_c}{(S_{i+1} + S_c)^2} + \alpha k {S_{i+1}}^{\alpha-1}\right) \end{aligned} \end{split}\]

Rubric:

  • 0.5 points if the derivative is correct

  • only 0.25 points if the derivative for one of the terms is incorrect, or if students started from the wrong expression

2.4 (0.5 point)

Provide the full expression used for iteration in the Newton-Raphson scheme, \(S_{i+1}^{j+1}\) (use your answer from task 2.3).

\[S_{i+1}^{j+1} = S_{i+1}^j - \frac{q(S_{i+1}^j)}{q^\prime(S_{i+1}^j)}\]
\[S_{i+1}^{j+1} = S_{i+1}^j - \frac{S_{i+1}^j - S_i - \Delta t \left(P(t_{i+1}) - \text{ET}_{\text{pot}}(t_{i+1}) \frac{S_{i+1}^j}{S_{i+1}^j + S_c} - k {S_{i+1}^j}^\alpha\right)}{1 + \Delta t \left(\text{ET}_{\text{pot}}(t_{i+1}) \frac{S_c}{(S_{i+1}^j + S_c)^2} + \alpha k {S_{i+1}^j}^{\alpha-1}\right)}\]

Rubric:

  • 0.5 points if the expression is fully correct

  • only 0.25 points if students started with the wrong expressions for \(q\) or \(q^\prime\) but correctly plugged them into the equation

Part 3#

3.1 (1 point)

Use a time step of 0.1 days, a tolerance of \(10^{-4}\) and a maximum number of Newton iteration of 50. What are the first five values of the time \(t\) and the storage \(S\) obtained with your implicit Euler scheme? Include at least four decimal digits.

time (days)

\(S\)

0.

10.

0.1

9.95680102

0.2

9.91373587

0.3

9.87080406

0.4

9.82800508

Rubric:

  • 1 point if all numbers are correct

  • only 0.5 points if the numbers for \(S\) are correct, but not correctly matched with the times.

3.2 (1 point)

Use a time step of 1 day, a tolerance of \(10^{-4}\) and a maximum number of Newton iteration of 50. Paste the plot from task 3.1 that compares your implicit Euler solution to the solution obtained with scipy.

How do the two solutions differ? How can you explain these differences?

  • The solution from the implicit Euler with relatively large time step is ahead in time of the Scipy solution.

  • This is because a purely implicit scheme uses the forcing data (rain, ET) of the next time step. If the time step is large, the scheme basically uses the data from one day ahead, explaining why the simulation runs ahead.

  • Note: In practice, this issue could be addressed by using an explicit scheme for the forcing term to get the timing right, but an implicit scheme for the remaining terms to make the scheme more stable.

Rubric:

  • 0.5 points for correctly describing the differences between the implicit scheme and the reference solution

  • additional 0.5 points for making the link between an implicit scheme and the mismatch in the timing of the forcing data

3.3 (1 point)

How does the tolerance relate to the number of Newton iterations needed? What can you conclude from this on the efficiency of the Newton scheme? (Use your results from task 3.3.)

The Newton iterations generally converge very quickly. Decreasing the tolerance by four orders of magnitude on average only requires about one additional iteration (about 2 iteration with \(\epsilon=10^{-2}\), 2.8 iterations with \(\epsilon=10^{-6}\) and 3 iterations with \(\epsilon=10^{-10}\)). This shows that the Newton scheme is quite efficient.

Rubric:

  • 0.25 points for more or less the right number of iterations (something between 1 and 3)

  • additional 0.5 points for noting an increase in the iterations with decreasing tolerance

  • additional 0.25 points for concluding that the Newton scheme is efficient

3.4 (2 points)

Reflect on the goodness of fit for the simulated discharge time series, including any temporal patterns in the deviations. What reasons could explain the misfit? What changes would you propose to improve the model?

How good is the model fit and how does this depend on time?

The model shows some patterns that can also be observed in the data:

  • precipitation events cause peaks in the discharge that last several days,

  • peaks in summer are lower than peaks in winter.

However, the overall model fit is poor. This is reasonable considering that the model is very simplistic. Some of the main shortcomings of the model are:

  • In summer, the discharge is generally too low, and sometimes almost drops to zero, whereas the observations show a baseflow in between peaks in summer.

  • The model overpredicts the discharge in the fall and early winter.

This pattern could indicate that the model is lacking a long-term storage term that takes up excess precipitation in the winter, stores it and discharges it as baseflow in the summer.

Potential reasons for the misfit and suggestions to improve the model

  • One main reason for a poor fit of the model to the data is that the representation of hydrological/physical processes in the model is too simplistic.

    • The model uses only a single storage term for the whole catchment. In reality, water is stored in many different compartments (vegetation, on the surface, unsaturated zone, groundwater, snow). Tom improve the model, one could add more storage components, and compose the total discharge of discharge from each storage. Differentiating a “slow” storage and a “fast” storage component could already help to better match seasonal patterns and baseflow.

    • The parametric equations used for the actual evapotranspiration and the storage-discharge relationship are empirical rather than based on physical principles. The relationships might not represent the actual relationships in this catchment well. Choosing a different parametric equation could potentially improve the model.

    • In the model, water can only leave the catchment through evapotranspiration or the river discharge. In reality, water could leave or enter the catchment through other points than the river, e.g. through groundwater. Then, the water balance that the model is based on is not correct.

  • Even if the model structure is appropriate, it might be that the model is not calibrated well. Adjusting the parameters might yield a better model fit.

  • The rain data used as input may not be representative for the catchment because it is from a single station. Using an average for the whole catchment (e.g., based on gridded data from a radar) instead of measurements from one point. In addition, the time resolution of the forcing data is only daily which might not be enough to resolve short term events.

Rubric:

  • 0.5 points if the answer provides a reasonable summary of the model fit (more descriptive than “good” or “poor”), e.g. mentioning specific aspects that look correct or not

  • additional 0.5 points if the answer discusses seasonal patterns in the model fit

  • additional 0.5 points if the answer mentions at least one reason for a poor model fit

  • additional 0.5 points if the answer makes at least one reasonable suggestion how to improve the model fit

  • up to 0.5 points can be subtracted for clearly wrong arguments

3.6 (1 point)

To get a more realistic representation of hydrologic processes, many models include several states. (For example, these states could represent a fast flow and slow flow component, or snow storage, but this is irrelevant for this question.) If you want to include such additional states in the model, how would this affect the numerical solution scheme and its computational cost (still using an implicit Euler scheme with Newton-Raphson iterations)?

When increasing the number of states, the unknown becomes a vector instead of a scalar. In the Newton-Raphson iterations…

  • …we then need to evaluate the Jacobian (the matrix of partial derivatives of the ODE with respect to the states) instead of a single derivative. For \(n\) states, this requires evaluating \(n^2\) derivatives.

  • …we need to invert the Jacobian in every iteration.

Overall, increasing the number of states increases the computational cost of the algorithm non-linearly. Solving an ODE with \(n\) states takes more than \(n\) times as long as solving an ODE with one state.

Rubric:

  • 1/3 point for simply mentioning that the computational effort increases.

  • additional 1/3 point for mentioning that we need to evaluate the Jacobian instead of just one derivative.

  • additional 1/3 point for discussing that the increase in required computational resources is more than linear with the number of states.

By Anna Störiko, Ronald Brinkgreve, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.

Contents