[BUG] FlexModel Calculated recharge problem for large gamma value

[HenrikSpa]

Calculating r using a large gamma causes to much reduction in the root zone storage state:

pastas/pastas/recharge.py

Lines 378 to 379 in 1cd9fd9

	# Calculate the recharge flux
	r[t] = min(ks * (sr[t] / srmax) ** gamma, sr[t])

The effect is that the root zone storage state gets drained more for large amount of rain and snowmelt than for smaller amounts. In other words, a small precipitation event causes the root zone to be a little wetter, while a large event drains it almost completely!

I made a semi-solution (there are probably better ways to do it) in the code below which divides the timestep into smaller iterative parts instead of one large step. The resulting image shows the problem with the regular calculation and the result using the iterative method.

An iterative step of 1 (the topmost plot) is equal to the regular calculation. One can see that when "sr before timestep" approaches and passes 150 the resulting "sr after timestep" starts to decrease for all gammas above 1.

For this example, using 5 steps almost solves the issue and there is very little, if any, difference between 5, 10 and 1000 iterative steps.

I'm not sure how often this problem exists and the model seems to work rather well anyway at the larger level and might be compensated in the response function, but it should probably be fixed somehow.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

mpl.use('Qt5Agg')

def regular(ks, gamma, srmax, sr_t):
    return min(ks * (sr_t / srmax) ** gamma, sr_t)

def iterative(ks, gamma, srmax, sr_t, steps):
    ks_per_step = ks / steps
    sr_sum = sr_t
    r = 0
    for step in range(steps):
        r_step = regular(ks_per_step, gamma, srmax, sr_sum)
        sr_sum -= r_step
        r += r_step
        if not sr_sum:
            break
    return r

ks = 150.0
srmax = 250.0
sr = [srmax*fill_level for fill_level in np.arange(0, 1.1, 0.1)]
test_number_of_iterative_steps = [1, 2, 5, 10, 1000]

fig, axs = plt.subplots(len(test_number_of_iterative_steps), 1, sharex='all', gridspec_kw={'hspace': 0.4,
                                                                                           'top': 0.95,
                                                                                           'right': 0.99,
                                                                                           'left': 0.1})
for ax_idx, steps in enumerate(test_number_of_iterative_steps):
    for gamma in [1, 4, 10, 20]:
        timesteps = list(range(len(sr)))
        sr_reg = [(sr[t] - regular(ks, gamma, srmax, sr[t])) for t in timesteps]
        sr_iter = [(sr[t] - iterative(ks, gamma, srmax, sr[t], steps=steps)) for t in timesteps]
        a, = axs[ax_idx].plot(sr, sr_reg, label=f'Regular gamma {gamma}', linestyle='-')
        axs[ax_idx].plot(sr, sr_iter, label=f'Iterative gamma {gamma}', linestyle='--', color=a.get_color())
        axs[ax_idx].set_title(f'{steps} iterative steps', fontsize=8)

axs[-1].legend(fontsize=6)
axs[-1].set_xlabel('sr before timestep')
fig.text(0.01, 0.5, 'sr after timestep', va='center', rotation='vertical')

plt.show()

Edit:
Here is an example from a calibration where this effect happened:

As can be seen in the image, the root zone storage is rather constant. During the first large snow melt the "Recharge (R)" is so great that it almost drains the root zone storage and during the large precipitation event it is completely drained.

[HenrikSpa]

But! Changing this behaviour will affect calibrated models that depend on the "flash drain"-effect of large amounts of precipitation and snowmelt.

[HenrikSpa]

My suggested solution doesn't completely solve even with above 1000 steps it seems.
EDIT: It does seems that 1000 solve the issue. (I made an error in the other code. The code in this issue seems correct).

Just so there is no misconception: When I though that the iterative solution didn't work, I was wrong. The iterative solution does seem to solve the original problem.

But as I wrote above, this change does have a large effect on the calibration. If this effect is used by the calibration to generate large amounts of recharge during heavy rains, the calibration will no longer work after the change. I'm now trying to offset this by adding a new parameter that controls whether excess precipitation is added to recharge (value 1.0) or removed as surface runoff (value 0.0).

[martinvonk]

We will look into the FlexModel and its parameters in #469. So duplicated.