diff --git a/docs/sphinx/source/whatsnew/v0.9.5.rst b/docs/sphinx/source/whatsnew/v0.9.5.rst
index 55ba6e3844..6c50f32642 100644
--- a/docs/sphinx/source/whatsnew/v0.9.5.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.5.rst
@@ -16,6 +16,14 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* Added the optional `string_factor` parameter to
+  :py:func:`pvlib.snow.loss_townsend` (:issue:`1636`, :pull:`1653`)
+
+Bug fixes
+~~~~~~~~~
+* Added a limit to :py:func:`pvlib.snow.loss_townsend` to guard against
+  incorrect loss results for systems that are near the ground (:issue:`1636`,
+  :pull:`1653`)
 
 * Added optional ``n_ar`` parameter to :py:func:`pvlib.iam.physical` to
   support an anti-reflective coating. (:issue:`1501`, :pull:`1616`)
@@ -53,6 +61,7 @@ Contributors
 ~~~~~~~~~~~~
 * Kevin Anderson (:ghuser:`kanderso-nrel`)
 * Will Holmgren (:ghuser:`wholmgren`)
+* Cliff Hansen (:ghuser:`cwhanse`)
 * Adam R. Jensen (:ghuser:`adamrjensen`)
 * Pratham Chauhan (:ghuser:`ooprathamm`)
 * Karel De Brabandere (:ghuser:`kdebrab`)
diff --git a/pvlib/snow.py b/pvlib/snow.py
index c06b317a2f..8cd1309bb6 100644
--- a/pvlib/snow.py
+++ b/pvlib/snow.py
@@ -219,7 +219,7 @@ def _townsend_effective_snow(snow_total, snow_events):
 
 def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
                   temp_air, poa_global, slant_height, lower_edge_height,
-                  angle_of_repose=40):
+                  string_factor=1.0, angle_of_repose=40):
     '''
     Calculates monthly snow loss based on the Townsend monthly snow loss
     model [1]_.
@@ -230,7 +230,8 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         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]_. [-]
+        Number of snowfall events each month. May be int or float type for
+        the average events in a typical month. Referred to as N in [1]_.
 
     surface_tilt : float
         Tilt angle of the array. [deg]
@@ -250,6 +251,11 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     lower_edge_height : float
         Distance from array lower edge to the ground. [m]
 
+    string_factor : float, default 1.0
+        Multiplier applied to monthly loss fraction. Use 1.0 if the DC array
+        has only one string of modules in the slant direction, use 0.75
+        otherwise. [-]
+
     angle_of_repose : float, default 40
         Piled snow angle, assumed to stabilize at 40°, the midpoint of
         25°-55° avalanching slope angles. [deg]
@@ -263,7 +269,12 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
     -----
     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``.
+    of ``surface_tilt``. The author of [1]_ recommends using one-half the
+    table width for ``slant_height``, i.e., the distance from the tracker
+    axis to the module edge.
+
+    The parameter `string_factor` is an enhancement added to the model after
+    publication of [1]_ per private communication with the model's author.
 
     References
     ----------
@@ -273,13 +284,22 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
        :doi:`10.1109/PVSC.2011.6186627`
     '''
 
+    # unit conversions from cm and m to in, from C to K, and from % to fraction
+    # doing this early to facilitate comparison of this code with [1]
+    snow_total_inches = snow_total / 2.54  # to inches
+    relative_humidity_fraction = relative_humidity / 100.
+    poa_global_kWh = poa_global / 1000.
+    slant_height_inches = slant_height * 39.37
+    lower_edge_height_inches = lower_edge_height * 39.37
+    temp_air_kelvin = temp_air + 273.15
+
     C1 = 5.7e04
     C2 = 0.51
 
-    snow_total_prev = np.roll(snow_total, 1)
+    snow_total_prev = np.roll(snow_total_inches, 1)
     snow_events_prev = np.roll(snow_events, 1)
 
-    effective_snow = _townsend_effective_snow(snow_total, snow_events)
+    effective_snow = _townsend_effective_snow(snow_total_inches, snow_events)
     effective_snow_prev = _townsend_effective_snow(
         snow_total_prev,
         snow_events_prev
@@ -288,37 +308,38 @@ def loss_townsend(snow_total, snow_events, surface_tilt, relative_humidity,
         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)
+    # the lower limit of 0.1 in^2 is per private communication with the model's
+    # author. CWH 1/30/2023
+    lower_edge_distance = np.clip(
+        lower_edge_height_inches**2 - effective_snow_weighted**2, a_min=0.1,
+        a_max=None)
     gamma = (
-        slant_height
-        * effective_snow_weighted_m
+        slant_height_inches
+        * effective_snow_weighted
         * cosd(surface_tilt)
-        / (lower_edge_height_clipped**2 - effective_snow_weighted_m**2)
+        / lower_edge_distance
         * 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
+        * effective_snow_weighted
         * cosd(surface_tilt)**2
         * ground_interference_term
         * relative_humidity_fraction
         / temp_air_kelvin**2
         / poa_global_kWh**0.67
+        * string_factor
     )
 
     return np.clip(loss_fraction, 0, 1)
diff --git a/pvlib/tests/test_snow.py b/pvlib/tests/test_snow.py
index 11f0b21e83..26ede9a782 100644
--- a/pvlib/tests/test_snow.py
+++ b/pvlib/tests/test_snow.py
@@ -6,6 +6,8 @@
 from pvlib import snow
 from pvlib.tools import sind
 
+import pytest
+
 
 def test_fully_covered_nrel():
     dt = pd.date_range(start="2019-1-1 12:00:00", end="2019-1-1 18:00:00",
@@ -108,6 +110,7 @@ def test__townsend_effective_snow():
 
 
 def test_loss_townsend():
+    # hand-calculated solution
     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])
@@ -118,6 +121,7 @@ def test_loss_townsend():
     poa_global = np.array([350000, 350000, 350000, 350000, 350000, 350000,
                            350000, 350000, 350000, 350000, 350000, 350000])
     angle_of_repose = 40
+    string_factor = 1.0
     slant_height = 2.54
     lower_edge_height = 0.254
     expected = np.array([0.07696253, 0.07992262, 0.06216201, 0.01715392, 0, 0,
@@ -125,5 +129,84 @@ def test_loss_townsend():
     actual = snow.loss_townsend(snow_total, snow_events, surface_tilt,
                                 relative_humidity, temp_air,
                                 poa_global, slant_height,
-                                lower_edge_height, angle_of_repose)
+                                lower_edge_height, string_factor,
+                                angle_of_repose)
     np.testing.assert_allclose(expected, actual, rtol=1e-05)
+
+
+@pytest.mark.parametrize(
+    'poa_global,surface_tilt,slant_height,lower_edge_height,string_factor,expected',  # noQA: E501
+    [
+        (np.asarray(
+            [60., 80., 100., 125., 175., 225., 225., 210., 175., 125., 90.,
+             60.], dtype=float) * 1000.,
+         2.,
+         79. / 39.37,
+         3. / 39.37,
+         1.0,
+         np.asarray(
+            [44, 34, 20, 9, 3, 1, 0, 0, 0, 2, 6, 25], dtype=float)
+         ),
+        (np.asarray(
+            [60., 80., 100., 125., 175., 225., 225., 210., 175., 125., 90.,
+             60.], dtype=float) * 1000.,
+         5.,
+         316 / 39.37,
+         120. / 39.37,
+         0.75,
+         np.asarray(
+            [22, 16, 9, 4, 1, 0, 0, 0, 0, 1, 2, 12], dtype=float)
+         ),
+        (np.asarray(
+            [60., 80., 100., 125., 175., 225., 225., 210., 175., 125., 90.,
+             60.], dtype=float) * 1000.,
+         23.,
+         158 / 39.27,
+         12 / 39.37,
+         0.75,
+         np.asarray(
+            [28, 21, 13, 6, 2, 0, 0, 0, 0, 1, 4, 16], dtype=float)
+         ),
+        (np.asarray(
+            [80., 100., 125., 150., 225., 300., 300., 275., 225., 150., 115.,
+             80.], dtype=float) * 1000.,
+         52.,
+         39.5 / 39.37,
+         34. / 39.37,
+         0.75,
+         np.asarray(
+             [7, 5, 3, 1, 0, 0, 0, 0, 0, 0, 1, 4], dtype=float)
+         ),
+        (np.asarray(
+            [80., 100., 125., 150., 225., 300., 300., 275., 225., 150., 115.,
+             80.], dtype=float) * 1000.,
+         60.,
+         39.5 / 39.37,
+         25. / 39.37,
+         1.,
+         np.asarray(
+             [7, 5, 3, 1, 0, 0, 0, 0, 0, 0, 1, 3], dtype=float)
+         )
+    ]
+)
+def test_loss_townsend_cases(poa_global, surface_tilt, slant_height,
+                             lower_edge_height, string_factor, expected):
+    # test cases from Townsend, 1/27/2023, addeed by cwh
+    # snow_total in inches, convert to cm for pvlib
+    snow_total = np.asarray(
+        [20, 15, 10, 4, 1.5, 0, 0, 0, 0, 1.5, 4, 15], dtype=float) * 2.54
+    # snow events are an average for each month
+    snow_events = np.asarray(
+        [5, 4.2, 2.8, 1.3, 0.8, 0, 0, 0, 0, 0.5, 1.5, 4.5], dtype=float)
+    # air temperature in C
+    temp_air = np.asarray(
+        [-6., -2., 1., 4., 7., 10., 13., 16., 14., 12., 7., -3.], dtype=float)
+    # relative humidity in %
+    relative_humidity = np.asarray(
+        [78., 80., 75., 65., 60., 55., 55., 55., 50., 55., 60., 70.],
+        dtype=float)
+    actual = snow.loss_townsend(
+        snow_total, snow_events, surface_tilt, relative_humidity, temp_air,
+        poa_global, slant_height, lower_edge_height, string_factor)
+    actual = np.around(actual * 100)
+    assert np.allclose(expected, actual)