diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index a69a3e18d4..f190b268e2 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -159,21 +159,25 @@ def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10):
     # swap angles so that phi[0,:,:,:] is the lesser angle
     phi.sort(axis=0)
 
-    # now re-use phi's memory again, this time storing cos(phi).
-    next_edge = phi[1, :, :, 1:]
-    np.cos(next_edge, out=next_edge)
-    prev_edge = phi[0, :, :, :-1]
-    np.cos(prev_edge, out=prev_edge)
-    # right edge of next row - left edge of previous row, again
-    # reusing memory so that the difference is stored in next_edge.
-    # Note that the 0.5 view factor coefficient is applied after summing
-    # as a minor speed optimization.
-    np.subtract(next_edge, prev_edge, out=next_edge)
-    np.clip(next_edge, a_min=0., a_max=None, out=next_edge)
-    vf = np.sum(next_edge, axis=-1) / 2
+    # compute cos(phi) into new arrays (avoid overwriting phi)
+    upper = phi[1].copy()
+    lower = phi[0].copy()
+    np.cos(upper, out=upper)
+    np.cos(lower, out=lower)
+
+    # include horizon bounds
+    # leftmost bound: cos(0)=1
+    # rightmost bound: cos(pi)=-1
+    upper = np.concatenate([upper, -np.ones_like(upper[..., :1])], axis=-1)
+    lower = np.concatenate([np.ones_like(lower[..., :1]), lower], axis=-1)
+
+    # wedge contribution = upper - lower for each wedge
+    diff = upper - lower
+    np.clip(diff, a_min=0., a_max=None, out=diff)
+
+    vf = np.sum(diff, axis=-1) / 2
     return vf
 
-
 def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10,
                            npoints=100, vectorize=False):
     """
diff --git a/tests/bifacial/test_utils.py b/tests/bifacial/test_utils.py
index e24af42b3e..e6344f1014 100644
--- a/tests/bifacial/test_utils.py
+++ b/tests/bifacial/test_utils.py
@@ -7,6 +7,7 @@
 from pvlib.shading import masking_angle, ground_angle
 from pvlib.tools import cosd
 from scipy.integrate import trapezoid
+from pvlib.bifacial.utils import vf_ground_sky_2d
 
 
 @pytest.fixture
@@ -190,3 +191,11 @@ def test_vf_ground_2d_integ(test_system_fixed_tilt):
     y = 0.5 * (1 - cosd(phi_y - ts['surface_tilt']))
     y1 = trapezoid(y, x) / (1 - 0)
     assert np.allclose(vf, y1, rtol=1e-2)
+
+def test_vf_ground_sky_2d_returns_valid_range():
+    vf = vf_ground_sky_2d(rotation=0, gcr=0.5, x=0.5,
+                          pitch=1, height=1, max_rows=3)
+
+    assert vf.shape == (1, 1)
+    assert np.all(vf >= 0.0)
+    assert np.all(vf <= 1.0)