Bug in PrairieInfil Implementation of Gray's 1985 frozen soil infiltration routine

[djmallum]

Summary

I believe there is some unintended behaviour in the control flow for the frozen soil part of the PrairieInfil module. The behaviour itself is not necessarily incorrect but its almost certainly not intended.

In short: If the first major melt occurs at the final melt day, no partitioning of the melt is done between rainfall or runoff. This could have unintended consequences. In addition, the code itself is a bit unclear.

Details

In this issue, I'll be referring to the model as "Crack".

PrarieInfil switches from Ayers for unfrozen soils to Crack for frozen soils when the SWE exceeds 25 mm. The module then waits for the first major melt, and then the Crack model is properly initialized by computing an infiltration index. The conditions for this initialization occur as follows:

1. No other major melts this season: crackstat[hh] = 0
2. It is a major melt: snowmelt[hh] >= Major[hh]
3. SWE is greater than 0. The specific check is SWE[hh] > Xinfil[2][hh] with the full code snippet:

crhmcode/crhmcode/src/modules/ClassPrairieInfil.cpp

Lines 216 to 224 in 61d577b

	else if (fallstat[hh] < 100.0)
	{
	if (snowmelt[hh] >= Major[hh] || crackstat[hh] >= 1)
	{
	if (SWE[hh] > Xinfil[2][hh] && snowmelt[hh] >= Major[hh])
	{
	infil_index(fallstat[hh] / 100.0, SWE[hh], Xinfil[0][hh], Xinfil[1][hh], infDays[hh]);
	Xinfil[2][hh] = SWE[hh];
	}

The Xinfil array (which is in dire need of renaming.... that aside) contains three values. The infiltration index as element 0, the maximum infiltrarion per major melt in element 1, and the SWE when Crack is first initialized at element 2. Since we are talking about pre-initialization, these values are all zero, as set here:

crhmcode/crhmcode/src/modules/ClassPrairieInfil.cpp

Lines 178 to 187 in 61d577b

	//If snow is larger than 25 mm and the frozen soil routine is off then we need to turn it on.
	if (SWE[hh] > 25.0 && !crackon[hh])
	{
	crackstat[hh] = 0;
	crackon[hh] = true; // initialise for current year
	timer[hh] = 0;
	Xinfil[0][hh] = 0.0;
	Xinfil[1][hh] = 0.0;
	Xinfil[2][hh] = 0.0;
	}

Therefore, there exists a hypothetical case where crack starts (i.e., SWE[hh] > 25.0 && !crackon[hh] returns true), then what proceeds is many minor melts (snowmelt[hh] >= Major[hh] returns false). Followed by a single large melt, eliminating the remaining snow cover. In this hypothetical case, the final melt will be registered as major correctly, however, all the values stored in Xinfil will be 0.

Infiltration is set according to:

crhmcode/crhmcode/src/modules/ClassPrairieInfil.cpp

Line 237 in 61d577b

snowinfil[hh] = snowmelt[hh] * Xinfil[0][hh];

, but as noted above, Xinfil[0][hh] is yet unset. Therefore, in this scenario, the amount of infiltrated melt for this day will be 0 and will always be runoff. In addition, any rain that might occur will also be runoff, rather than infiltrated. And during the melt period, rain is certainly possible.

Right now, due to the large threshold for "frozen soils" at 25 mm, this is unlikely to occur. However, if the threshold were to be lowered or allowed to be set by the user in future versions, this bug could create some anomalous runoff at the end of the melt season (which is short enough already).

Expected behaviour

Even if the first major melt is the final melt, the index is still evaluated and melt is partitioned. A smarter system would have to account for these minor melts. Defaulting to runoff is dangerous, especially because the soil module already will divert excess infiltrated amounts already if the soil is saturated.

Suggested change

The quickest and easiest change is simple to change line 220 from

crhmcode/crhmcode/src/modules/ClassPrairieInfil.cpp

Line 220 in 61d577b

if (SWE[hh] > Xinfil[2][hh] && snowmelt[hh] >= Major[hh])

to

if (SWE[hh] >= Xinfil[2][hh] && snowmelt[hh] >= Major[hh])

The >= operation ensures that the index is computed in the unique case of SWE = 0 and the melt is partitioned between runoff and infiltration. At the start of the melt period, when Xinfil[2]=0.0, SWE is guaranteed to be positive. This should only impact the very last melt in a season with no major melts. Even though this is potentially rare, behaviour like this being implicit (and not directly handled) is a bad idea in terms of code clarity.