Error in snow.loss_townsend with respect to gamma quantity?

[smmeredith]

I'm testing the snow.loss_townsend function, and I noticed it was giving me different results than what I expected.

After looking through the code, it seems the major reason is due to how gamma is calculated.

gamma = (
        slant_height_inches
        * effective_snow_weighted
        * cosd(surface_tilt)
        / lower_edge_distance
        * 2
        * tand(angle_of_repose)
    )

In the pvlib function, the quantity tand(angle_of_repose) is in the numerator of the quotient. In the Townsend paper, the formula is given as gamma = [R*Se’cos(tilt)]/[(H2 – Se’2)/2tan(P)], which makes the same quantity a denominator.

It's a little ambiguous since the equation is poorly formatted in the original publication I'm referencing, but still seems to imply that this tangent should be in the denominator.

Can the function author confirm that this is an error?

[cwhanse]

As published that equation is equivalent to

$$ \frac{R Se' \cos(tilt)}{\frac{H^2 - Se'^2}{2 \tan(P)}} $$

so the pvlib code is what was published. I agree; the bracket placement can be a bit challenging. There's a derivation of the denominator (mass of snow dissipated) in this comment, so I think the code represents the paper author's intent.

[cwhanse]

Can you elaborate on the results you expected to see?

[smmeredith]

Thanks for the reply Cliff.

I'm comparing results generated from my own Python code which calls the snow.loss_townsend function to results from an Excel file which uses the Townsend snow loss model (not created by me, BTW, which is partly why I'm doing this comparison). I'm trying to understand the differences between the Excel model and the algorithm in the function. There are a few differences, but most significant seems to be how gamma is calculated.

Regarding the comment that you linked to, I'm actually not super familiar with the Townsend model, but if the denominator of gamma is supposed to be a summation of those areas, then that derivation makes sense.

My problem is with how the equation for gamma is formatted in the original publication. The equation for gamma (eq. 5) appears there as:

If there were parentheses around the 2*tan(P) term, that would be consistent with the Python code. But the lack of parentheses makes it seem like tan(P) should overall be in the denominator of the gamma quotient.

[kandersolar]

@smmeredith Having looked into the model in some depth while making that diagram Cliff linked, I am convinced that pvlib's choice of whether to multiply or divide by that tangent is consistent with the model's intent. So if the doubt here stems from the fact that the gamma equation's formatting leaves room for interpretation, I suggest thinking about this "ratio of snow received/snow dissipated" quantity, chewing on the diagram and accompanying math, and seeing if you become similarly convinced of which side of the fraction that tangent ought to be on.

But if there is a discrepancy between pvlib's implementation and some more authoritative implementation of the model, that would certainly be interesting. Any details you can provide regarding whatever calculation difference you are seeing would be appreciated!

[smmeredith]

@kandersolar
To be clear, I think I agree on the pvlib implementation as well (given the derivation of the formula) - I was just commenting on why the formula in the paper was ambiguous.

Regarding the differences between the pvlib function and the Excel file, it was really a number of things - different formula for gamma, having a minimum number of snow events to calculate effective snowfall, different weighting for effective snowfall, , etc. But after looking at the Townsend publication, the pvlib function, and the Excel implementation of the model, I think I agree with the pvlib code where there are any discrepancies.

[mikofski]

Thanks for following up ❤️. It’s good to raise concerns and question implementations.

[allison-mueller]

I came here looking for the same thing - thanks @smmeredith for confirming I'm not crazy. The formula does imply that the tan() term is in the denominator!

[kandersolar]

Perhaps a note on this topic in the function's docstring is warranted. Pull requests are welcome!