Add some COSERO functionalities to HydPy-H-Land

[tyralla]

This is a follow-up to #67. At the time of writing, the first version of our HBV96-PREVAH combination seems to work well. Before performing detailed comparisons, we prefer to add an analogous application model that combines process equations of HBV96 and COSERO. As for PREVAH, we take the modules "above the soil" from HBV96 and the modules "below the soil" from COSERO.

[tyralla]

COSERO offers different recession constants for base flow (TAB3) and lake outflow (TABLA). Hence, surface water bodies and groundwater seem to be more disconnected than in HBV96. Nevertheless, to keep the modelling of water bodies consistent among our HBV derivatives, I again start with considering them as HBV96-like "internal lakes" that are directly connected with the lower groundwater storage.

[tyralla]

Both the surface runoff and the interflow module should behave highly non-linear due to the applied thresholds. Therefore, I start to calculate them on the zone level and follow the PREVAH instead of HBV96 discretisation.

This decision is consistent with the table given in section 11.2.2 of the project report "Spatio-temporal water balance Danube". The same table says that COSERO calculates the additional runoff concentration in the "subbasin route reservoir" also on the zone level. However, if its implementet linearly, this should not be be necessary. Still an open question...

[tyralla]

Besides the mentioned thresholds, all runoff concentration modules seem to follow the linear storage approach. Principally, we could also allow for additional nonlinearity by implementing a nonlinearity parameter like HBV96's Alpha. But then, we would not be able to apply the approximate analytical solutions described in COSERO's documentation. Also, calibration would possibly become too complicated. Hence, I also start with implementing them linearly. Also, it seems reasonable to use the given analytical solutions.

[tyralla]

The implementation of the analytical solution of the two-outlet reservoirs is tricky and requires much testing. The following code allows verifying that the simple explicit Euler method reproduces the results documented in the docstrings of method Calc_QAb1_QVs1_BW1_V1 and method Calc_QAb2_QVs2_BW2_V1 (with much more computational effort and less accuracy):

from numpy import inf
steps = 500000
digits = 4
dt = 1.0 / steps
nul = 0.0
for test in range(10):
    h = [4.0, 4.0, 4.0, 4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0][test]
    k1 = [2.0, inf, nul, 2.0, 2.0, 2.0, inf, nul, 2.0, 2.0][test]
    k2 = [4.0, 4.0, 4.0, inf, nul, 4.0, 4.0, 4.0, inf, nul][test]
    for hru in range(8):
        s0 = [2.0, 2.0, 9.0, 9.0, 2.0, 5.0, 2.0, 6.0][hru]
        qz = [0.0, 2.0, 0.0, 2.0, 5.0, 0.1, 0.5, 2.5][hru]
        if test == 0:
            qa1_test = round([0.0, 0.0, 1.561119, 1.956437, 0.246114, 0.181758, 0.0, 1.0][hru], digits)
            qa2_test = round([0.442398, 0.672805, 1.780559, 1.978219, 1.007586, 1.086805, 0.5, 1.5][hru], digits)
            s0_test = round([1.557602, 3.327195, 5.658322, 7.065344, 5.7463, 3.831437, 2.0, 6.0][hru], digits)
        if test == 1:
            qa1_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            qa2_test = round([0.442398, 0.672805, 1.990793, 2.221199, 1.018414, 1.117516, 0.5, 1.615203][hru], digits)
            s0_test = round([1.557602, 3.327195, 7.009207, 8.778801, 5.981586, 3.982484, 2.0, 6.884797][hru], digits)
        if test == 2:
            qa1_test = round([0.0, 0.0, 5.0, 6.0, 2.115471, 1.0, 0.0, 3.5][hru], digits)
            qa2_test = round([0.442398, 0.672805, 0.884797, 1.0, 0.884529, 0.896317, 0.5, 1.0][hru], digits)
            s0_test = round([1.557602, 3.327195, 3.115203, 4.0, 4.0, 3.203683, 2.0, 4.0][hru], digits)
        if test == 3:
            qa1_test = round([0.0, 0.0, 1.967347, 2.393469, 0.408182, 0.414775, 0.0, 1.319592][hru], digits)
            qa2_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            s0_test = round([2.0, 4.0, 7.032653, 8.606531, 6.591818, 4.685225, 2.5, 7.180408][hru], digits)
        if test == 4:
            qa1_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            qa2_test = round([2.0, 4.0, 9.0, 11.0, 7.0, 5.1, 2.5, 8.5][hru], digits)
            s0_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
        if test == 5:
            qa1_test = round([0.703511, 1.09883, 3.165801, 3.561119, 1.691807, 1.778544, 0.802341, 2.604682][hru], digits)
            qa2_test = round([0.351756, 0.549415, 1.5829, 1.780559, 0.845904, 0.889272, 0.40117, 1.302341][hru], digits)
            s0_test = round([0.944733, 2.351756, 4.251299, 5.658322, 4.462289, 2.432184, 1.296489, 4.592977][hru], digits)
        if test == 6:
            qa1_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            qa2_test = round([0.442398, 0.672805, 1.990793, 2.221199, 1.018414, 1.117516, 0.5, 1.615203][hru], digits)
            s0_test = round([1.557602, 3.327195, 7.009207, 8.778801, 5.981586, 3.982484, 2.0, 6.884797][hru], digits)
        if test == 7:
            qa1_test = round([2.0, 4.0, 9.0, 11.0, 7.0, 5.1, 2.5, 8.5][hru], digits)
            qa2_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            s0_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
        if test == 8:
            qa1_test = round([0.786939, 1.213061, 3.541224, 3.967347, 1.852245, 1.988653, 0.893469, 2.893469][hru], digits)
            qa2_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            s0_test = round([1.213061, 2.786939, 5.458776, 7.032653, 5.147755, 3.111347, 1.606531, 5.606531][hru], digits)
        if test == 9:
            qa1_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
            qa2_test = round([2.0, 4.0, 9.0, 11.0, 7.0, 5.1, 2.5, 8.5][hru], digits)
            s0_test = round([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0][hru], digits)
        s0_ = s0
        qa1, qa2 = 0.0, 0.0
        for _ in range(steps):
            s0 += dt * qz
            if s0 > h:
                qa1_ = s0 - h
                if k1 > 0.0:
                    qa1_ *= min(dt / k1, 1.0)
            else:
                qa1_ = 0.0
            qa2_ = s0
            if k2 > 0.0:
                qa2_ *= min(dt / k2, 1.0)
            qa1 += qa1_
            qa2 += qa2_
            s0 -= qa1_ + qa2_
        bal = round(s0_ + qz - qa1 - qa2 - s0, digits)
        assert bal == 0.0, f"test={test}, hru={hru}: {bal} != 0.0"
        qa1 = round(qa1, digits)
        assert qa1 == qa1_test, f"test={test}, hru={hru}: {qa1} != {qa1_test}"
        qa2 = round(qa2, digits)
        assert qa2 == qa2_test, f"test={test}, hru={hru}: {qa2} != {qa2_test}"
        s0 = round(s0, digits)
        assert s0 == s0_test, f"test={test}, hru={hru}: {s0} != {s0_test}"

.

[tyralla]

COSERO models the baseflow reservoir as a linear runoff concentration storage and uses the analytical solution for constant inflow to solve it. HBV96 (at least our L-Land Version 1) relies on the same equation but uses the explicit Euler method for solving it. Baseflow reacts slowly, so switching between different integration methods should not make a big difference here. Hence, for now, I decide to reuse the related method of L-Land Version 1 for our HBV96/COSERO implementation. At least, this helps to avoid duplicate code. However, we should discuss the names of the different state and flux sequences of all of our HBV-like models later to increase their consistency.

[tyralla]

Due to the decision above: The original COSERO recession coefficients TAb1 and TAb2 (as well as TVs1 and TVs2) have the unit "T", while the original HBV96 recession coefficient K4 has the unit "1/T". This inconsistency can complicate working with the model, for example when trying to set restrictions like TAb1 < TAb2 < 1/K4.