Snow submodels

@BGWKlein

We are going to tackle the next big step in rolling out the submodel concept by moving the snow modules from the existing land models into the HydPy-Snow model family.

There are likely many things to discuss. This issue will serve as a notebook to keep track of our decisions.

This change could affect the following main models:

gland_gr4 (case GR)
gland_gr5 (case GR)
gland_gr6 (case GR)
hland_96 (case HBV)
hland_96c (case HBV)
hland_96p (case HBV)
lland_dd (case DD+)
lland_knauf (case Knauf)
lland_knauf_ic (case Knauf & Ic)
whmod_rural (case DD)
whmod_urban (case DD)
wland_wag (case WALRUS)
wland_gd (case WALRUS)

Unfortunately, we have to deal with six different "cases":

GR: These main models possess no own snow module so far. Hence, we only need to add the functionality of using snow submodels. The other unique aspect is that all GR models are lumped, but that should pose no problem.
HBV: All HBV models work with the same snow module so far. Basically, this utilises the simple day-degree method (including minor extensions such as refreezing the snow layer's liquid water content). However, there are also structurally more relevant features, such as reflecting local spatial heterogeneity by defining snow classes and redistributing snow from higher to lower zones to prevent “snow towers” in high-mountain areas.
DD+: The day-degree method without snow classes and a different "snow tower" prevention mechanism.
Knauf: Without snow classes and without a "snow tower" prevention mechanism, but with much more complex energy balance calculations. This includes an energy exchange between the snow layer and the underlying soil, which might further complicate matters.
Knauf & Ic: Even more complex, because snow on the ground and snow in the canopies are modelled separately.
DD: Very similar to the DD+ case, but even simpler, because there is no separation between frozen and liquid water in the snow layer (and thus no refreezing) and also no "snow tower" prevention mechanism.
WALRUS: Conceptually identical to the DD case. However, the HydPy-W models are the only main models defined by ODE systems rather than approximate ad hoc solutions to the original equations.

The first stack of questions:

Do we want to support the WALRUS case? (It would be cool, of course. But can we, for example, formulate the Knauf method as an ODE system with a reasonable amount of work?)
Do we want to support the Knauf & Ic case? (We already started to discuss this point. Our tendency was not to support it, and maybe even to remove lland_knauf_ic, because we are not overly happy with its functionality. The documentation states it is "still under discussion".)
Should there be only one snow submodel interface? Or do we want to separate, for example, snow melt and snow retention? Similar thought for "local" snow processes and snow redistribution.

Also, we must be aware of causing a major breaking change (except for the GR case). I am afraid it would be too much work to duplicate all affected models during a transition phase (e.g. add hland_96_new now, and remove hland_96 and rename hland_96_new to hland_96 later). But maybe we can find some tricks that simplify converting existing projects. (Would it be possible to add default submodels that work like the old snow modules, and to direct values passed to the old parameters to them?)

cc @BGWKlein

Do we want the snow submodels to be optional? Implementing this feature should not be much additional work. On the other hand, I don't see a clear benefit to ignoring snow processes. (At least the simple day-degree method should never cause much computational overhead.)

To make things more concrete (and to identify further points for discussion), I will begin by introducing the simplest possible snow submodel, which I will call HydPy-Snow-DD (Degree-day-based snow modelling) for now, and making it usable by whmod_rural.

Should the snow submodels also be applicable as standalone models? This could be a more interesting feature than making them optional and may simplify testing.

@sivogel

Spatial resolution

HydPy-GR is lumped. All other land models are semi-distributed.

Most snow modules distinguish only between different hydrological response units. The snow module of HydPy-H further separates all HRUs into the same number of snow classes. The snow module of HydPy-Snow-CN divides a lumped subcatchment (like one modelled with HydPy-GR) into layers (different heights).

To avoid creating too many submodels, I will add the concept of snow classes to HydPy-Snow-DD. Setting the number of snow classes to 1 leads to the simple "DD case".

Regarding the layer approach, I am unsure. I see the following options:

Make sure it is always only used in combination with HydPy-GR or with a single HRU.
Allow that the main model is semi-distributed, but the snow submodel is still lumped (except for the layering). This would require some aggregation (of precipitation) and disaggregation (of water release across different layers).
Allow for a layering of the individual HRUs.

Option 3 would have the benefit of enforcing the same dimensionality, as the layer-dimension replaces the snow class-dimension. On the downside, it would possibly make little sense if the HRUs themselves are based on elevation classes (up to the user). Also, the configuration could become complex (not sure, I have little experience with HydPy-Snow so far).

cc @sivogel

I would prefer Option 3 because it provides the greatest flexibility. E.g. using option H-ZONEN (elevation zones) in LARSIM, the snow pack of each Hydrological Response Unit (HRU) is divided into the same number of elevation layers based on the subbasin's elevation distribution. While the area fractions of the different snow layers can vary from one layer to the next, the configuration, meaning the number of layers, the mean elevation, and the area fraction of each individual layer, is identical across all HRUs.

Implementing Option 3 would perfectly allow HydPy to accommodate this LARSIM-like snow layer approach (specifically matching the Option SCHNEE: H-ZONEN EXTERN configuration).
This option could also be used for the HBV setup, in this case the area fractions and mean elevations would simply be identical.

To implement this efficiently, I recommend enforcing the following structural restrictions across all HRUs:

The number of layers must be identical for each HRU.
The area fraction of each corresponding layer must be identical for each HRU (with the fractions of the layers summing up to 1.0), though individual layer fractions can differ from one another.
The mean elevation of each layer must be identical for each HRU.

Yes, I agree that the flexibility of appealing. However, maybe I understand the proposed restrictions incorrectly, but wouldn't "identical for each HRU" result in the following limitations?

We could not reflect the absence of woods on high mountains. (This would require, e.g., different layerings for field and forest zones.)
We could not totally follow the HBV setup (at least as we implemented it in HydPy-H), because here the mean elevation of each layer (snow class) is identical to the HRU's elevation.
For HRUs with the same landcover type, the input from the snow submodel would often be very similar or even identical.

We could not reflect the absence of woods on high mountains. (This would require, e.g., different layerings for field and forest zones.)

Yes, this is exactly how it is implemented in LARSIM, and it is indeed a inherent structural limitation of this option. Depending on the spatial scale and heterogeneity of your subbasins, this approach can have some drawbacks. However, allowing layers to be defined individually for each separate HRU introduces too much configuration complexity and overhead. In cases where this vegetation-elevation mismatch becomes a significant issue, the cleanest hydrological workaround is simply to define smaller, more homogeneous subbasins.

We could not totally follow the HBV setup (...) because here the mean elevation of each layer (snow class) is identical to the HRU's elevation.

Yes, you are completely right. This specific case needs to be explicitly accounted for in the model logic.

To bridge this gap without breaking the structural uniformity of Option 3, we can implement a fallback rule: if no explicit layer elevations are defined in the configuration, the model should automatically default to using the HRU's own elevation for all layers.

We had an offline discussion on the above matter and concluded that we should first focus on the more fundamental question of the scope of individual model types. Put simply: Who is responsible for reading air temperature measurements (or other meteorological factors) and eventually adjusting them to different elevations?

After thinking it over, it appears to me that we need to move more functionality into the submodels. Currently, for example, HydPy-H (and other land models) reads the subbasins' average air temperature and adjusts it to individual heights. It then uses the adjusted values to calculate snowmelt and passes them to the HydPy-Evap submodels. After the submodularisation of snow processes, HydPy-H no longer requires temperature data for its own equations.

I see the following benefits for removing the temperature-handling from HydPy-H (and similar land models) and moving it to HydPy-Evap and HydPy-Snow (or, eventually, some sub-submodels):

HydPy-H would not need to read input data that it does not require for itself.
We would get rid of the inconsistency that HydPy-H provides its submodels with temperature data but not, e.g., with wind data (which depends also on elevation and is required by some evaporation and snow process equations).
We would gain some flexibility regarding modelling evaporation and snow processes. For example, we could calculate snowmelt using adjusted temperature data, but calculate reference evapotranspiration based on subbasin temperature data and adjust the resulting values afterwards.
We would gain some flexibility regarding modelling evaporation or snow with different approaches. For example, we could implement a snow submodel that adjusts temperature to the elevation of individual zones (HBV style) or to subbasin-wide layers (GR style and eventually LARSIM style).

We already discussed the last point and would prefer not to create submodels that differ in both this spatial aspect and their actual snowmelt equations (to avoid too many combinations). Maybe there are ways to prevent this by separating into "HydPy-Snow-Spatial-Structure" submodels and "HydPy-Snow-Process-Equation" sub-submodels. Or by creating a new type of submodel for adjusting (input) data to individual hydrological response units (that could be freely used by HydPy-Evap and HydPy-Snow).

A protection mechanism to prevent accidental connections. (Maybe we could add an ALLOWS_CROSSCONNECTIONS flag or something more specific to the submodels.)

This can be done later and is likely related to the required consistency checks for submodel cross-connection. I opened #194 as a reminder.

Water areas

More methods of whmod_rural than expected relied on special assumptions about how the implemented snow routines function. A special case is the snowmelt threshold temperature, which was represented by the literal value 0.0, which controlled the "ponding".

The handling of water areas is a similar issue but relevant to (nearly) all land models. Currently, no implemented snow routine accounts for the possibility of snow covering a frozen lake. The land models know this and sometimes redirect the water flow around the snow routine. I suggest we let all land models pass the throughfall to and query the release from the new snow submodels for all land types. Then, it is up to the snow submodel to decide whether to simulate snow processes for water areas or simply return the received throughfall.

Receiving input

So far, I have misused the Get_Precipitation_V1 function in whmod_rural by letting it return simulated throughfall rather than measured precipitation. It requires a little work, but it seems consistent to also add the submodel interfaces ThroughfallModel_V1 (which defines only get_throughfall - for using main models as sub-submodels) and ThroughfallModel_V2 (which defines some more methods - for real submodels), to let all land models inherit from ThroughfallModel_V1, and to implement an IO submodel that follows ThroughfallModel_V2.

1-layer submodels would only call get_throughfall to get their input (snowpack on the ground). 2-layer submodels would also call get_precipitation (snowpack in the canopy).

The above decision regarding water areas has the consequence that the whmod flux sequence Ponding now also gets a snow submodel's water release for water areas, despite being described as "Ponding on land surfaces". "Ponding on land or water areas" sounds strange. Maybe we should rename the sequence. In lland, it's WaDa, and in hland, it's In_. (It seems no one can come up with a good name for precipitation, throughfall or snowmelt receiving the land or water surface.)

Glaciers

The combination of whmod_rural and whmod_urban with our preliminary version of snow_dd (which so far only has the capabilities required for WHMod) works satisfactorily (no changes in the integration test results). Now I am starting to extract the snow routines from our HydPy-H model family. I hope there are no further complicated technical problems, but we will need to discuss some hydrological aspects. One such thing is the handling of glaciers, which are implemented as special land types both in HydPy-H and HydPy-L.

Do we want to

leave the glacier-related routines in HydPy-H (and HydPy-L),
move the glacier-related routines to HydPy-Snow,
create special submodels in the new model family HydPy-Glacier to be used by the HydPy-H (and HydPy-L) main models,
or create special sub-submodels in the new model family HydPy-Glacier to be used by the HydPy-Snow submodels?

We need to take into account that the definition of a glacier zone can trigger the partial conversion of snow into glacier ice. If this conversion is to be calculated by the glacier routines, a glacier submodel would need to modify the states of a snow submodel. This aspect speaks against option 3, because it would extend the scope of submodel cross-connections beyond simple exchanges of state information.

At first glance, option 4 could be the most favourable for flexibility and simplicity. It also has the benefit that we can start with option 1 and can postpone introducing HydPy-Glacier to a later step.

However, if we want to include the glacier routines in HydPy-Snow, we should do so right away to avoid unnecessary work (the inter-model communication between the snow and glacier modules).

Temperature-specific precipitation correction

An interesting aspect of all HydPy-H application models is that they allow for correcting precipitation with respect to the relative amounts of liquid and frozen water (which makes sense because wind-caused underestimations of precipitation are usually more severe for snowflakes than for raindrops). This correction occurs before precipitation hits the interception storage. Hence, I will leave the parameters TT and TTInt in hland and implement their analogues in snow_dd. This duplication is a little unfortunate. Maybe we could improve the situation by passing the values automatically when adding the snow submodel to the hland main model.

(Isolating and moving this part of the whole, relatively complex precipitation adjustment of HydPy-H to HydPy-Snow does not seem right to me. Especially since I consider it likely that we will also create precipitation-adjustment submodels at some point.)

The temperature-specific precipitation correction of HydPy-H is strange. In my words, FracRain represents the fraction of initially guessed snowfall and measured total precipitation. If we use a higher factor for adjusting snowfall than for rainfall, the fraction of the corrected snowfall and the corrected total precipitation becomes smaller than FracRain (assuming FracRain is not zero or one). This means, if we let HydPy-H pass the throughfall (which might be identical to corrected total precipitation) to HydPy-Snow, it cannot calculate the liquid and frozen parts of throughfall by estimating FracRain as HydPy-H did, because FracRain would then represent the "real" fraction of snowfall and total precipitation.

I see the following options:

Give HydPy-Snow the necessary knowledge about the HydPy-H internals. (Bad idea.)
Let HydPy-H pass the liquid and frozen parts of throughfall instead of their sum. (Also no good idea - I think the snow submodels should do this division.)
Just don't care that FracRain means slightly different things in HydPy-H and HydPy-Snow and accept some changes in the results. (Maybe only very few people would recognise this inconsistency, but if one starts thinking about the details, it's confusing.)
Adjust HydPy-H's calculation of FracRain so that it represents the "real" fraction of snowfall and total precipitation. (I did not check if this could work, but I think so; in my view, this would increase consistency and understandability, but would also cause changes in simulation results, of course.)
Just remove the whole temperature-specific precipitation correction from HydPy-H. (I never set the rain- and the snow-related coefficients differently - this could be handled in pre-processing; however, it means we would remove a feature).

We decided to implement option 4.

We can change the equation

$$ P_c = Rain_c + Snow_c = Rain \cdot RfCF + Snow \cdot SfCF = FracRain \cdot P \cdot RfCF + (1 - FracRain) \cdot P \cdot SfCF $$

to

$$ P = Rain + Snow = \frac{FracRain \cdot P_c}{RfCF} + \frac{(1 - FracRain) \cdot P_c}{SfCF} $$

which can be rearranged to

$$ P_c = \frac{P}{\frac{FracRain}{RfCF} + \frac{1 - FracRain}{SfCF}} $$

For example, setting $FracRain=0.2$, $RfCF=1.1$, and $SfCF=1.3$ would result in a total (weighted) precipitation correction of $1.2543859649122806$ instead of $1.26$.

For example, setting $FracRain=0.2$, $RfCF=1.1$, and $SfCF=1.3$ would result in a total (weighted) precipitation correction of $1.2543859649122806$ instead of $1.26$.

Under some circumstances, the effect of the planned change is larger than I expected. Taking up the last example, the amounts of rain and snow would be $0.2 \cdot 1.1 = 0.22$ and $0.8 \cdot 1.3 = 1.04$ mm in the current state and $0.2 \cdot 1.2543859649122806 = 0.25087719298245614$ and $0.2 \cdot 1.2543859649122806 = 1.0035087719298246$. Without an initial snow layer, without any melting, and with a typical water-holding capacity ($WHC$) of 0.1, the release ($In$) is $0.22 - 0.1 \cdot 1.04 = 0.116$ and $0.25087719298245614 - 0.1 \cdot 1.0035087719298246 = 0.15052631578947367$, respectively. This is a difference of 30 %.

Adjust the temperature-sensitive precipitation correction of all appl…

We decided to accept the described difference. I added the following note to the documentation of method Calc_PC_V1:

Apply the precipitation correction factors and adjust precipitation to the altitude of the individual zones.

Basic equation:

$PC = \frac{P \cdot PCorr \cdot \big(1 + PCAlt \cdot (ZoneZ - Z) \big)} {\frac{FracRain}{RfCF} + \frac{1 - FracRain}{SfCF}}$

Until HydPy 6, Calc_PC_V1 relied on the following basic equation:

$PC = P \cdot PCorr \cdot (1 + PCAlt \cdot (ZoneZ - Z)) \cdot (RfC + SfC)$

The fraction of $RfC$ and $RfC + SfC$ was often not consistent with $FracRain$. $FracRain$ did represent the relative amount of estimated, uncorrected rainfall with respect to measured total precipitation. Since HydPy 7, $FracRain$ always represents the relative amount of estimated, corrected rainfall to corrected total precipitation. This might not be totally consistent with HBV96 but makes reasoning and the coupling of snow submodels simpler. See issue 188 for background information.

I implemented the change in commit 05d033f. This commit also contains many changes (mostly tiny ones) in some doctests that explain central HydPy concepts, because these often rely on the Lahn example project, which is based on hland_96.

Glaciers

Do we want to

leave the glacier-related routines in HydPy-H (and HydPy-L),
move the glacier-related routines to HydPy-Snow,
create special submodels in the new model family HydPy-Glacier to be used by the HydPy-H (and HydPy-L) main models,
or create special sub-submodels in the new model family HydPy-Glacier to be used by the HydPy-Snow submodels?

We decided on the first option (at least temporarily). I will thus add parameter Cov to HydPy-H (SC is already taken) because the existence of a snow cover prevents glacier melting. (And I will likely add a similar parameter to HydPy-L later.)

Update `hydpy.core.modeltools.Model.find_submodels`.

snow_dd now also works in combination with hland_96, hland_96p, and hland_96c. This makes it necessary to adjust our Lahn example project and much of our documentation, including many doctests. I expect this will reveal some additional problems we need to handle. The first one that I encountered is that the parameter values of snow_dd are now written twice in each control file, and its state values are now written twice in each condition file - one time due to its role as "real submodel" of hland_96 and the other time due to evap_aet_hbv96 building a cross connection to it to query its snow cover. The underlying problem is that find_submodels returns snow_dd for both of these roles.

I suggest adding a parameter named include_crossconnections which defaults to False. find_submodels should usually return snow_dd for those submodel ports where it has been "really" added (via an add_snowsubmodel_v1 statement). With this convention, automatically written control and condition files look plausible, and a single snow_dd is not visited twice when we, for example, want to read its time series.

However, I am not sure if there are some places where we would want to visit snow_dd twice (by setting Include_crossconnections to True). So far, the only (not overly important) place I found is the string representation of hland_96.Model. Sticking to the default, it might look like this:

    hland_96
        aetmodel: evap_aet_hbv96
            petmodel: evap_pet_hbv96
        rconcmodel: rconc_uh
        snowmodel: snow_dd

After setting Include_crossconnections to True in hydpy.core.modeltools.Model.__repr__, the result might look so:

    hland_96
        aetmodel: evap_aet_hbv96
            petmodel: evap_pet_hbv96
            snowcovermodel: snow_dd
        rconcmodel: rconc_uh
        snowmodel: snow_dd

Personally, I prefer the second option, but I am open to suggestions. Maybe we should add an explanation like this:

    hland_96
        aetmodel: evap_aet_hbv96
            petmodel: evap_pet_hbv96
            snowcovermodel: snow_dd (cross connection)
        rconcmodel: rconc_uh
        snowmodel: snow_dd

Navigation Menu

Repository navigation

Snow submodels #188

Spatial resolution

6 remaining items

Water areas

Receiving input

Glaciers

Temperature-specific precipitation correction

Glaciers

Update `hydpy.core.modeltools.Model.find_submodels`.

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Snow submodels #188

Description

Activity

tyralla commented on Jun 2, 2026

tyralla commented on Jun 2, 2026

tyralla commented on Jun 2, 2026

tyralla commented on Jun 2, 2026

Spatial resolution

BGWKlein commented on Jun 2, 2026

tyralla commented on Jun 2, 2026

BGWKlein commented on Jun 3, 2026

tyralla commented on Jun 9, 2026

6 remaining items

tyralla commented on Jul 9, 2026

tyralla commented on Jul 10, 2026

Water areas

tyralla commented on Jul 10, 2026

Receiving input

tyralla commented on Jul 10, 2026

tyralla commented on Jul 15, 2026

Glaciers

tyralla commented on Jul 15, 2026

Temperature-specific precipitation correction

tyralla commented on Jul 18, 2026

tyralla commented on Jul 20, 2026

tyralla commented on Jul 21, 2026

tyralla commented on Jul 21, 2026

tyralla commented on Jul 22, 2026

Glaciers

tyralla commented on Sep 1, 2026

Update hydpy.core.modeltools.Model.find_submodels.

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Update `hydpy.core.modeltools.Model.find_submodels`.