Townsend snow (building on #1251)

docs/sphinx/source/reference/effects_on_pv_system_output.rst

@@ -21,6 +21,7 @@ Snow
    snow.coverage_nrel
    snow.fully_covered_nrel
    snow.dc_loss_nrel
+   snow.loss_townsend
 
 Soiling
 -------

docs/sphinx/source/whatsnew/v0.9.2.rst

@@ -18,6 +18,8 @@ Bug fixes
   where passing localized timezones with large UTC offsets could return
   rise/set/transit times for the wrong day in recent versions of ``ephem``
   (:issue:`1449`, :pull:`1448`)
+* Added Townsend-Powers monthly snow loss model: :py:func:`pvlib.snow.loss_townsend`
+  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 
 
 Testing
@@ -40,3 +42,5 @@ Contributors
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
 * Chencheng Luo (:ghuser:`roger-lcc`)
 * Prajwal Borkar (:ghuser:`PrajwalBorkar`) 
+* Abhishek Parikh (:ghuser:`abhisheksparikh`)
+* Taos Transue (:ghuser:`reepoi`)

pvlib/snow.py

@@ -5,7 +5,7 @@
 
 import numpy as np
 import pandas as pd
-from pvlib.tools import sind
+from pvlib.tools import sind, cosd, tand
 
 
 def _time_delta_in_hours(times):
@@ -185,3 +185,130 @@ def dc_loss_nrel(snow_coverage, num_strings):
        Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
     '''
     return np.ceil(snow_coverage * num_strings) / num_strings
+
+
+def _townsend_effective_snow(snow_load, snow_events):
+    '''
+    Calculates effective snow using the total snowfall in inches received
+    each month and the number of snowfall events each month.
+
+    Parameters
+    ----------
+    snow_load : array-like
+        Inches of snow received each month. Referred to as S in the paper [in]
+
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in the paper [-]
+
+    Returns
+    -------
+    effective_snowfall : array-like
+        Effective snowfall as defined in the Townsend model
+
+    References
+    ----------
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
+       update from two winters of measurements in the SIERRA. Conference
+       Record of the IEEE Photovoltaic Specialists Conference.
+       003231-003236. :doi:`10.1109/PVSC.2011.6186627`
+       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
+    '''  # noqa: E501
+    snow_events_no_zeros = np.maximum(snow_events, 1)
+    effective_snow = 0.5 * snow_load * (1 + 1 / snow_events_no_zeros)
+    return np.where(snow_events > 0, effective_snow, 0)
+
+
+def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
+                  temp_air, poa_global, slant_height, lower_edge_drop_height,
+                  angle_of_repose=40):
+    '''
+    Calculates monthly snow loss based on the Townsend monthly snow loss model [1]_.
+
+    Parameters
+    ----------
+    snow_total : array-like
+        Inches of snow received each month. Referred to as S in the paper [in]
+
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in the paper [-]
+
+    surface_tilt : float
+        Tilt angle of the array [deg]
+
+    relative_humidity : array-like
+        Monthly average relative humidity [%]
+
+    temp_air : array-like
+        Monthly average ambient temperature [C]
+
+    poa_global : array-like
+        Monthly plane of array insolation [kWh/m2/month]
+
+    slant_height : float
+        Row length in the slanted plane of array dimension [in]
+
+    lower_edge_drop_height : float
+        Drop height from array edge to ground [in]
+
+    angle_of_repose : float, default 40
+        piled snow angle, assumed to stabilize at 40°, the midpoint of
+        25°-55° avalanching slope angles [deg]
+
+    Returns
+    -------
+    loss : array-like
+        Monthly average DC capacity loss fraction due to snow coverage
+
+    Notes
+    -----
+    This model has not been validated for tracking arrays; however, for tracking
+    arrays [1]_ suggests using the maximum rotation angle in place of surface_tilt.
+
+    References
+    ----------
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
+       update from two winters of measurements in the SIERRA. Conference
+       Record of the IEEE Photovoltaic Specialists Conference.
+       003231-003236. 10.1109/PVSC.2011.6186627.
+       Available at https://www.researchgate.net/publication/261042016_Photovoltaics_and_snow_An_update_from_two_winters_of_measurements_in_the_SIERRA
+    '''  # noqa: E501
+
+    C1 = 5.7e04
+    C2 = 0.51
+
+    snow_total_prev = np.roll(snow_total, 1)
+    snow_events_prev = np.roll(snow_events, 1)
+
+    effective_snow = _townsend_effective_snow(snow_total, snow_events)
+    effective_snow_prev = _townsend_effective_snow(
+        snow_total_prev,
+        snow_events_prev
+    )
+    effective_snow_weighted = (
+        1 / 3 * effective_snow_prev
+        + 2 / 3 * effective_snow
+    )
+
+    drop_height_clipped = np.maximum(lower_edge_drop_height, 0.01)
+    gamma = (
+        slant_height
+        * effective_snow_weighted
+        * cosd(surface_tilt)
+        / (drop_height_clipped**2 - effective_snow_weighted**2)
+        * 2
+        * tand(angle_of_repose)
+    )
+
+    ground_interference_term = 1 - C2 * np.exp(-gamma)
+    temp_air_kelvin = temp_air + 273.15
+    loss_percentage = (
+        C1
+        * effective_snow_weighted
+        * cosd(surface_tilt)**2
+        * ground_interference_term
+        * relative_humidity
+        / temp_air_kelvin**2
+        / poa_global**0.67
+    )
+
+    return loss_percentage / 100

pvlib/tests/test_snow.py

@@ -95,3 +95,32 @@ def test_dc_loss_nrel():
     expected = pd.Series([1, 1, .5, .625, .25, .5, 0])
     actual = snow.dc_loss_nrel(snow_coverage, num_strings)
     assert_series_equal(expected, actual)
+
+
+def test__townsend_effective_snow():
+    snow_load = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    expected = np.array([7.5, 7.5, 5, 0, 0, 0, 0, 0, 0, 0, 3.75, 6.66666667])
+    actual = snow._townsend_effective_snow(snow_load, snow_events)
+    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+
+
+def test_loss_townsend():
+    snow_total = np.array([10, 10, 5, 1, 0, 0, 0, 0, 0, 0, 5, 10])
+    snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    surface_tilt = 20
+    relative_humidity = np.array([80, 80, 80, 80, 80, 80, 80, 80, 80, 80,
+                                  80, 80])
+    temp_air = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+    poa_global = np.array([350, 350, 350, 350, 350, 350, 350, 350, 350, 350,
+                           350, 350])
+    angle_of_repose = 40
+    slant_height = 100
+    lower_edge_drop_height = 10
+    expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
+                         0, 0, 0, 0, 0.02643821, 0.06068194])
+    actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
+                                relative_humidity, temp_air,
+                                poa_global, slant_height,
+                                lower_edge_drop_height, angle_of_repose)
+    np.testing.assert_allclose(expected, actual, rtol=1e-05)