Finish the implementation of GR4J and friends

[tyralla]

@BGWKlein implemented prototypes for GR4J and some related models a long time ago, and @sivogel did some polishing a year ago in #112 and recently updated her PR to account for the new submodel concept (regarding evaporation.

It might be a little late, but I open this issue as a notebook to keep track of the final discussions.

[tyralla]

• Rename grxj_land. I propose to rename the model group grxjland to gr_land (gr for genie rurale).

This suggestion is taken from the PR. There seems to be no fundamental difference between the model versions for daily and hourly simulation steps. Hence, committing the "j" is reasonable. Then, the "x" is not required as a placeholder anymore. However, we only took the first letters from HBV96, LARSIM, and WALRUS (hland, lland, and wland) to make clearer distinctions between the original and the HydPy implementations. So, the alternative suggestion is gland. However, note that a suggested name for one of the evap models is currently evap_pet_hbv96. What do others think?

[BGWKlein]

I would prefer gr_land as model group name.

…

-----Ursprüngliche Nachricht----- Von: Christoph Tyralla ***@***.***> Gesendet: Dienstag, 16. April 2024 09:39 An: hydpy-dev/hydpy ***@***.***> Cc: Klein, Bastian, M2, B40 ***@***.***>; Mention ***@***.***> Betreff: Re: [hydpy-dev/hydpy] Finish the implementation of GR4J and friends (Issue #134) * Rename grxj_land. I propose to rename the model group grxjland to gr_land (gr for genie rurale). This suggestion is taken from the PR <#112 (comment)> . There seems to be no fundamental difference between the model versions for daily and hourly simulation steps. Hence, committing the "j" is reasonable. Then, the "x" is not required as a placeholder anymore. However, we only took the first letters from HBV96, LARSIM, and WALRUS (hland, lland, and wland) to make clearer distinctions between the original and the HydPy implementations. So, the alternative suggestion is gland. However, note that a suggested name for one of the evap models is currently evap_pet_hbv96. What do others think? — Reply to this email directly, view it on GitHub <#134 (comment)> , or unsubscribe <https://github.com/notifications/unsubscribe-auth/AOZHFRAEQ55RNJ7DO7EGKLLY5TISXAVCNFSM6AAAAABGIY2EBKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDANJYGQZTQMRYGU> . You are receiving this because you were mentioned. <https://github.com/notifications/beacon/AOZHFRA7K26YJBMCCGOUNTDY5TISXA5CNFSM6AAAAABGIY2EBKWGG33NNVSW45C7OR4XAZNMJFZXG5LFINXW23LFNZ2KUY3PNVWWK3TUL5UWJTT2WFDI2.gif> Message ID: ***@***.***>

[tyralla]

@sivogel also raised the question (in an offline discussion) of whether to extract the unit hydrograph into a submodel (see #130). Currently, @parvumflumen implements a suitable submodel family in #133. I see no reason against this (except about one or two days of additional work).

Please note the related discussion in #130.

[tyralla]

I would prefer gr_land as model group name.

@sivogel: When (or if) you apply this name, could you please carefully check if the underscore does not cause any problems? Code sections might assume an application model's family name to be the first substring before the underscore. I would most likely expect such assumptions in the code that generates the documentation.

[tyralla]

The current proposal suggests three application models with limited differences. Can we reduce the number of application models without losing functionality? This could help remove some redundancies (especially in testing). However, we should refrain from such refactoring if it affects comprehensibility (too much).

[tyralla]

I would prefer gr_land as model group name.

I remembered one place where this name could cause problems: autodoc_applicationmodel (see the split operation; I guess we explain the "conventional way" somewhere else)

[tyralla]

In an offline discussion, we decided on gland.

[sivogel]

Currently we have 3 different gland application models gland_gr4, gland_gr5 and gland_gr6. @tyralla asked if it is possible to merge these 3 application models into one application model. In my opinion we should stay with 3 application models because I do not see a way to merge Calc_F_V1 and Calc_F_V2 into 1 application model without having meaningless parameters.

• we keep 3 application models for gland

But maybe we can combine the snow model with hysteresis effect with the one without. To do this, we would need to introduce GLocalMax into the model without hysteresis effect. This would be an overhead in the simpler application model, but would save us two application models.

• we combine those 2 application models

[sivogel]

The method Calc_PLayer_V1 calculates the precipitation for each height level of the snow module. In airGR it is assumed that above 4000m precipitation is no longer height dependent. This assumption was not found in any other publication. Should we keep the airGR assumption?

• We decided to stick with the 4000 m precipiation threshold

[sivogel]

I could not find any differences between the hourly and the daily gr4j airGR model. The Fortran scripts on airGR differ only in the parameter NH, which specifies the length of the unit hydrograph outflow. However, this is later truncated to the parameter x4 (unless x4 is greater than 20 days (480h), in which case it is truncated to 20 days (480h)). Therefore, in my current opinion, no separate implementation of GR4H is necessary.

• we found more differences:
• The exponent in the daily unit hydrograph is 2.5 in daily model and 1.25 in hourly model (we add the exponent as optional parameter to the uh submodel and set it by default to 2.5)
• in gr5h an interception model is introduced (we try to make the interception module optional to all models)

[sivogel]

The actual evaporation is calculated from the potential evaporation using equation: Calc_Es_V1 based on the filling of the production store. In contrast to the previous evaporation approaches, the net evaporation (i.e. pet - p) is used in the calculation instead of the potential evaporation, for which a suitable submodel interface would first have to be created that also passes on the precipitation or only the net evaporation to the submodel. This should not be too difficult to implement, but requires a little programming effort

• the calculation of the actual evapotranspiration of gland is very specific, we decided not to submodularize

[sivogel]

gr5j in the air package airGR splits the output of the unit hydrograph, but in some publications gr5j (like gr4j and gr6j) the input to the unit hydrograph is split and two seperate unit hydrographs are applied. Are we sure that we want to stay with the airGR version?
Lavenne:

airGR:

• We stick with the airGR version

[sivogel]

For equation Calc_SolidFraction_V1 the literature Turcotte2007 is referenced, but no page is given. The document has over 400 pages and no specific page is given so we cannot prove the correctness of this method.
For equation Calc_SolidFraction_V1 the threshold of 1500m does not seem to make sense because it leads to a unsteady fuction (see doctest). The difference in the results for 1500m depends strongly on the distribution of Tmin, Tmax and Tmean. It is therefore not possible to find a continuous transition.

• Bastian Klein will ask airGR developers about this

[sivogel]

Currently all variable names are taken from the airGR package. Should we adjust them. Should we find more descriptive names or leave then as they are.

• we stay with the old ones

[sivogel]

The equation
$$Melt = ((1 - MinMelt) \cdot GRatio + MinMelt) \cdot PotMelt$$
leads to small melt rates for small amounts of snow. In summer the amount of snow is therefore approaching 0, but never reaches 0. Might this be a problem for us and should we fix this?

• there should not be any conflicts, but we decided to set the snow amount to zero if there is nearly no snow. Therefore we introduce a threshold value

[sivogel]

We assume (similarly to the implementation of hland) that zinputs is not the station height, but always the mean height of the catchment area. Therefore the zinputs parameter is no longer necessary and can be calculated by zlayers and layerarea. We decided to remove the hypsodata option from zlayers and add a pre-processing method to the model to calculate layerarea and zlayers from hypsodata, as layerarea should not be set independently of hypsodata.

[sivogel]

I have updated my comments with the results of our discussions yesterday and today.

[sivogel]

In the current HydPy-gland version we have a mixture of parameter and series names from:

• airGR R Code
• airGR Fortran Code
• airGR documentation (graphics)

In the following I want to give an overview which can help us to

HydPy (current)	Fortran	R	documentation	HydPy (proposed)
Area
X1	Param(1)	x1	X1	X1
X2	Param(2)	x2	X2	X2
X3	Param(3)	x3	X3	X3
X4	Param(4)	x4	X4	X4
X5	Param(5)	x5	X5	X5
X6	Param(6)	x6	X6	X6
IMax	Imax	Imax	Imax	IMax
P	P1/Precip	Precip	P	P
PET	E/PE	PotEvap	E	E
En	EN		En	EN
Pn	PN/Pn	Pn	Pn	PN
Ps	PS/Ps	Ps	Ps	PS
EI	EI	EI	Ei	EI
Es	ES		Es	ES
AE	AE	AE		AE
Pr	PR	PR	Pr	PR
PrUH1	PRHU1			PR9
PrUH2	PRHU2			PR1
QOutUH2	StUH2			Q10
Perc	PERC/Perc	Perc	Perc	Perc
Q9	StUH1(1)/Q9	Q9	Q9	Q9
Q1	StUH2(1)/Q1	Q1	Q1	Q1
	EXCH/Exch	Exch	F	FPot
FD	AEXCH2/Aech2	Exch2	F	FD
FR	AEXCH1/Aexch1	Exch1	F	FR
	AEXCH1+AEXCH2/AExch	AExch		FTot
FR2	EXCH	Exch	F	FR2
Qr	QR	QR	Qr	QR
Qr2	QRExp	QRExp	QrExp	QR2
Qd	QD	QD	Qd	QD
Qt	Q/Qsim	Qsim	Q	QH
				QV
I	St(3)	Interc	I	I
S	St(1)/Prod	Prod	S	S
R	St(2)/Rout	Rout	R	R
R2	St(3)	Exp	R2	R2

Feel free to comment and make suggestions

[BGWKlein]

Shouldn't we use for input Sequences like Precipitation the standard HydPy names? I can't find the issue mentioning these names. Is there a standard name for potential evaporation?
I can't find the variable mean temperature in the list.

[tyralla]

We discussed this in #75. We now support multiple naming conventions. The actual implementation of an InputSequence subclass (its name) usually follows the model-specific nomenclature. The HydPy standard name is added via the class variable STANDARD_NAME.

The current list of the standard names is here. We added no "Mean_" prefixes to the names (so far) because reading averages is the standard case while, for example, reading maximums is special.

[tyralla]

The potential evapotranspiration standard name is defined here (I overlooked it at first sight. We should better sort all names in strict alphabetical order.)

[sivogel]

We still have problems in making our model available in different time steps. Ficchi 2019 (doi:10.1016/j.jhydrol.2019.05.084 ) describes in Table A.1 that the equations for

• percolation (perc)
• discharge from recharge (qr)
• groundwater exchange (f)
  must be adapted for different time steps. Le Moine shows in his PhD thesis 2008 (Appendix C) that this is due to different results of analytical integration for different time steps.
  However, in the airGR package only the percolation equation is adapted for different time steps. I suspect that this is only done for percolation, because for routing storage and groundwater exchange this adaptation can be "hidden" in the calibration parameters x2 and x3 (but would conflict with the given unit of x2 and x3 in the documentation). For percolation, a "hidden" adjustment of x1 is not possible, since x1 is also used in the calculation of ES and PS and was therefore included in the equations.
  Currently we do not fully understand the derivation of the analytical integration in Le Moines PHD thesis.

[tyralla]

The GR community seems to have put some effort into changing the simulation time step without changing the simulation result (too much). This even includes writing GR versions in a continuous state-space form (which would be relatively simple to implement in HydPy but is currently not on our target list).

I find it relatively hard to find all relevant information, partly because I am unfamiliar with French. However, as our current goal is "only" to achieve an implementation that is functionally similar to airGR, we probably do not need to understand all the details of the underlying assumptions. As far as I can see, all equations of these GR versions rely on the analytical integration of individual fluxes (operator splitting). Also, there is always the problem with the assumed temporal distribution of meteorological input, especially precipitation, when dealing with long simulation steps. So, we cannot expect that equivalent daily and hourly parameters result in equivalent daily and hourly simulation results (the same holds, e.g., for HydPy-H and HydPy-L). In my opinion, it is sufficient if our GR implementations are as good as those of arGR, except that we should ideally support arbitrary simulation step widths, not only daily and hourly ones.

I try to adjust the relevant equations of our current GR code accordingly and note my assumptions here.

[tyralla]

percolation

In Ficchi (2017), I found the following equation for calculating the percolation constant:

$$ \beta = 5.25 \cdot (3600 / \Delta t)^{1/4} $$

For $\Delta t = 1h$, this is precisely $21/4$ like in GRxH, but for $\Delta t = 1d$, it is only approximately $9/4$ like in GRxJ.

The following examples demonstrate that the above equation is correct (identical simulation results when applying the percolation equation on a daily and an hourly basis):

def perc(*, s: float, x1: float, v: float) -> float:
    return s * (1.0 - (1.0 + (v * s / x1) ** 4.0) ** (-0.25))


s0 = 10.0
x1 = 30.0

perc_1d = perc(s=s0, x1=x1, v=4.0 / 9.0)
assert round(perc_1d, 8) == 0.00120391, perc_1d

perc_1d = perc(s=s0, x1=x1, v=1.0 / (5.25 * (1.0 / 24.0) ** (1.0 / 4.0)))
assert round(perc_1d, 8) == 0.00097482, perc_1d

s = s0
perc_1h_cum = 0.0
for _ in range(24):
    perc_1h_now = perc(s=s, x1=x1, v=4.0 / 21.0)
    s -= perc_1h_now
    perc_1h_cum += perc_1h_now
assert round(perc_1h_cum, 8) == 0.00097482, perc_1h_cum

Unfortunately, we cannot apply the above equation without losing consistency with GRxJ. I am only speculating, but the ratio $9 / 4$ could be simply due to rounding, so I suggest adding the percolation constant as a "derived parameter" to our GR implementation and calculating its value automatically with the following in equation:

def calc_beta(dt: float) -> None:
    print(f"{round(4.0 * (5.25 * (1.0 / dt) ** 0.25))} / 4")

for hours in (0.25, 0.5, 1.0, 3.0, 6.0, 12.0, 24.0):
    calc_beta(dt=hours)
# 30 / 4
# 25 / 4
# 21 / 4
# 16 / 4
# 13 / 4
# 11 / 4
# 9 / 4

[tyralla]

groundwater exchange and outflow of the nonlinear routing store

Regarding the outflow term $QR$, things are in principle like discussed above, except that Ficchi suggests modifying the model parameter $x3$ directly instead of an inaccessible, time-dependent model constant. Maybe the only reason for this difference is that $x1$ is used within multiple equations and $x3$ only in the single equation which is affected by time step issues. It could also be that airGR does not implement these conversions because it just expects all parameter values to be given in agreement with the relevant time step (1 day or 1 hour). But I am only speculating. As HydPy separates between parameterstep and simulationstep and provides automatic conversions if possible, we should implement them, and I suggest doing so by overriding the get_timefactor method of the Parameter base class. The following example confirms (or at least does not contradict) that the conversion function works accurately:

def calc_qr(*, r: float, x3: float) -> float:
    return r * (1.0 - (1.0 + (r / x3) ** 4.0) ** -0.25)


def calc_fr(*, r: float, x2: float, x3: float) -> float:
    return x2 * (r / x3) ** 3.5


x2_1d = 7.0
x3_1d = 30.0
x2_1h = x2_1d * (1.0 / 24.0) ** 0.125
x3_1h = x3_1d * (1.0 / 24.0) ** -0.25

r0 = 10.0

qr_1d = calc_qr(r=r0, x3=x3_1d)
assert round(qr_1d, 8) == 0.03062823, qr_1d
r = r0
qr_1h_cum = 0.0
for _ in range(24):
    qr_1h_now = calc_qr(r=r, x3=x3_1h)
    r -= qr_1h_now
    qr_1h_cum += qr_1h_now
assert round(qr_1h_cum, 8) == 0.03062823, qr_1h_cum  # okay

fr_1d = calc_fr(r=r0, x2=x2_1d, x3=x3_1d)
assert round(fr_1d, 8) == 0.14968340, fr_1d
r = r0
fr_1h_cum = 0.0
for _ in range(24):
    fr_1h_now = calc_fr(r=r, x2=x2_1h, x3=x3_1h)
    r -= fr_1h_now
    fr_1h_cum += fr_1h_now
assert round(fr_1h_cum, 8) == 0.14602887, fr_1h_cum  # ???

For the groundwater exchange $FR$, there is a (in this case) tiny difference between the results of the daily and the hourly simulation. I cannot exclude a misunderstanding on my side, but I suspect the given transformation is only an approximation. After a short look at the relevant sections, neither Perrin nor Ficchi mentions the original differential equation, and Ficchi only gives some notes on the derivation of the other transformation equations.

Additionally, for GRx5 and GRx6, the calculation of $FR$ also includes parameter $x5$, which might result in further differences.

However, the given equation (and our current understanding) is the best we have, so I will implement the transformation for $x2$ as described for $x3$. If someone knows how to improve things, we can do it later.

[tyralla]

Regarding IMax, I encountered an interesting paragraph in the airGR docs:

Function which calculates the maximal capacity of the GR5H interception store. This function compares
the interception evapotranspiration from the GR5H interception store for different maximal
capacity values with the interception evapotranspiration classically used in the daily GR models
(e.g. GR4J). Among all the TestedValues, the value that gives the closest interception evapotranspiration
flux over the whole period is kept.

I am unsure if we ever need such functionality, so I will ignore it for now but post it as a reminder for potential later discussions.

[tyralla]

Reviewing all the GR methods (equations) again, I found two water balance mistakes. Update_R_V1 and Update_R_V2 use r = max(r + fluxes, 0) statements but update none of the fluxes in case r + fluxes would be negative. I checked this for gland_gr4, and check_waterbalance can, in fact, return non-zero values for negative values of X4.

I remember there was one place in the documentation where one exchange term was said to be "potential". I removed this word because it was not consistent with other descriptions. If FD, FR, and FR2 are used in the respective check_waterbalance methods (as they are), they must be "actual".

If FR is "actual" instead of "potential", we must reconsider the execution order of the methods of gland_gr6. Currently, Update_R_V2 (which also should update FR) is called before Calc_FR2_V1, which just defines fr2 = fr. We should change this order so as not to let our pure water balance-related fix affect the actual discharge simulation results of gland_gr6.

[tyralla]

While fixing the above problems, I noticed one more thing. It appears to me that when switching from GR5 to GR6, the routing storage-related groundwater exchange doubles:

$$ 2 \cdot FR_{GR5} = FR_{GR6} + FR2_{GR6} $$

And later it is, for example, $0.4 \cdot Q9 + FR2$ instead of $0.4 \cdot (Q9 + FR2)$.

I will leave it like this, but someone more familiar with GR6 may reconsider this.

[tyralla]

I will leave it like this, but someone more familiar with GR6 may reconsider this.

It may be okay because the routing store receives less input. So, it is usually less filled and $FR$ is, on average, smaller for GR6 than for GR4 and GR5.

[tyralla]

I created a new branch, removed all files related to HydPy-Snow, opened PR #160, and squash merged it into the master branch.

So, HydPy-G is finally ready and the online documentation should be updated soon. See https://hydpy-dev.github.io/hydpy/master/HydPy-GA.html.

This is the branch for finishing HydPy-Snow: https://github.com/hydpy-dev/hydpy/tree/feature/models/snow