diff --git a/docs/sphinx/source/reference/effects_on_pv_system_output.rst b/docs/sphinx/source/reference/effects_on_pv_system_output.rst
index 92efa946b4..62c1035860 100644
--- a/docs/sphinx/source/reference/effects_on_pv_system_output.rst
+++ b/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
 -------
diff --git a/docs/sphinx/source/whatsnew/v0.9.3.rst b/docs/sphinx/source/whatsnew/v0.9.3.rst
index 2f3c643d05..3bd6c45e4f 100644
--- a/docs/sphinx/source/whatsnew/v0.9.3.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.3.rst
@@ -12,6 +12,9 @@ Enhancements
 * Low resolution altitude lookup map
   :py:func:`~pvlib.location.lookup_altitude`
   (:issue:`1516`, :pull:`1518`)
+* Added Townsend-Powers monthly snow loss model:
+  :py:func:`pvlib.snow.loss_townsend`
+  (:issue:`1246`, :pull:`1251`, :pull:`1468`)
 
 Bug fixes
 ~~~~~~~~~
@@ -37,3 +40,5 @@ Contributors
 ~~~~~~~~~~~~
 * João Guilherme (:ghuser:`joaoguilhermeS`)
 * Nicolas Martinez (:ghuser:`nicomt`)
+* Abhishek Parikh (:ghuser:`abhisheksparikh`)
+* Taos Transue (:ghuser:`reepoi`)
diff --git a/pvlib/snow.py b/pvlib/snow.py
index bf40c7b995..c06b317a2f 100644
--- a/pvlib/snow.py
+++ b/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,140 @@ 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_total, snow_events):
+    '''
+    Calculates effective snow using the total snowfall received each month and
+    the number of snowfall events each month.
+
+    Parameters
+    ----------
+    snow_total : array-like
+        Snow received each month. Referred to as S in [1]_. [cm]
+
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in [1]_. [-]
+
+    Returns
+    -------
+    effective_snowfall : array-like
+        Effective snowfall as defined in the Townsend model. [cm]
+
+    References
+    ----------
+    .. [1] Townsend, Tim & Powers, Loren. (2011). Photovoltaics and snow: An
+       update from two winters of measurements in the SIERRA. 37th IEEE
+       Photovoltaic Specialists Conference, Seattle, WA, USA.
+       :doi:`10.1109/PVSC.2011.6186627`
+    '''
+    snow_events_no_zeros = np.maximum(snow_events, 1)
+    effective_snow = 0.5 * snow_total * (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_height,
+                  angle_of_repose=40):
+    '''
+    Calculates monthly snow loss based on the Townsend monthly snow loss
+    model [1]_.
+
+    Parameters
+    ----------
+    snow_total : array-like
+        Snow received each month. Referred to as S in [1]_. [cm]
+
+    snow_events : array-like
+        Number of snowfall events each month. Referred to as N in [1]_. [-]
+
+    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. [Wh/m2]
+
+    slant_height : float
+        Row length in the slanted plane of array dimension. [m]
+
+    lower_edge_height : float
+        Distance from array lower edge to the ground. [m]
+
+    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. 37th IEEE
+       Photovoltaic Specialists Conference, Seattle, WA, USA.
+       :doi:`10.1109/PVSC.2011.6186627`
+    '''
+
+    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
+    )
+    effective_snow_weighted_m = effective_snow_weighted / 100
+
+    lower_edge_height_clipped = np.maximum(lower_edge_height, 0.01)
+    gamma = (
+        slant_height
+        * effective_snow_weighted_m
+        * cosd(surface_tilt)
+        / (lower_edge_height_clipped**2 - effective_snow_weighted_m**2)
+        * 2
+        * tand(angle_of_repose)
+    )
+
+    ground_interference_term = 1 - C2 * np.exp(-gamma)
+    relative_humidity_fraction = relative_humidity / 100
+    temp_air_kelvin = temp_air + 273.15
+    effective_snow_weighted_in = effective_snow_weighted / 2.54
+    poa_global_kWh = poa_global / 1000
+
+    # Calculate Eqn. 3 in the reference.
+    # Although the reference says Eqn. 3 calculates percentage loss, the y-axis
+    # of Figure 7 indicates Eqn. 3 calculates fractional loss. Since the slope
+    # of the line in Figure 7 is the same as C1 in Eqn. 3, it is assumed that
+    # Eqn. 3 calculates fractional loss.
+    loss_fraction = (
+        C1
+        * effective_snow_weighted_in
+        * cosd(surface_tilt)**2
+        * ground_interference_term
+        * relative_humidity_fraction
+        / temp_air_kelvin**2
+        / poa_global_kWh**0.67
+    )
+
+    return np.clip(loss_fraction, 0, 1)
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 3db56c8f61..11f0b21e83 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -95,3 +95,35 @@ 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_total = np.array([25.4, 25.4, 12.7, 2.54, 0, 0, 0, 0, 0, 0, 12.7,
+                           25.4])
+    snow_events = np.array([2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 2, 3])
+    expected = np.array([19.05, 19.05, 12.7, 0, 0, 0, 0, 0, 0, 0, 9.525,
+                         254 / 15])
+    actual = snow._townsend_effective_snow(snow_total, snow_events)
+    np.testing.assert_allclose(expected, actual, rtol=1e-07)
+
+
+def test_loss_townsend():
+    snow_total = np.array([25.4, 25.4, 12.7, 2.54, 0, 0, 0, 0, 0, 0, 12.7,
+                           25.4])
+    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([350000, 350000, 350000, 350000, 350000, 350000,
+                           350000, 350000, 350000, 350000, 350000, 350000])
+    angle_of_repose = 40
+    slant_height = 2.54
+    lower_edge_height = 0.254
+    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_height, angle_of_repose)
+    np.testing.assert_allclose(expected, actual, rtol=1e-05)