diff --git a/docs/examples/shading/plot_passias_diffuse_shading.py b/docs/examples/shading/plot_passias_diffuse_shading.py
index 828cf141a1..751c1bcc40 100644
--- a/docs/examples/shading/plot_passias_diffuse_shading.py
+++ b/docs/examples/shading/plot_passias_diffuse_shading.py
@@ -20,15 +20,21 @@
 # :py:func:`pvlib.shading.masking_angle_passias` and
 # :py:func:`pvlib.shading.sky_diffuse_passias`.
 #
+# However, the pvlib-python authors believe that this approach is incorrect.
+# A correction is suggested and compared with the diffuse shading as obtained
+# with the view factor model.
+#
 # References
 # ----------
 #  .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
 #     panels", Solar Cells, Volume 11, Pages 281-291.  1984.
 #     DOI: 10.1016/0379-6787(84)90017-6
 
-from pvlib import shading, irradiance
 import matplotlib.pyplot as plt
 import numpy as np
+from cycler import cycler
+
+from pvlib import bifacial, shading, irradiance
 
 # %%
 # First we'll recreate Figure 4, showing how the average masking angle varies
@@ -43,7 +49,7 @@
 
 plt.figure()
 for k in [1, 1.5, 2, 2.5, 3, 4, 5, 7, 10]:
-    gcr = 1/k
+    gcr = 1 / k
     psi = shading.masking_angle_passias(surface_tilt, gcr)
     plt.plot(surface_tilt, psi, label=f'k={k}')
 
@@ -60,17 +66,18 @@
 # diffuse plane of array irradiance (after accounting for shading) to diffuse
 # horizontal irradiance. This means that the deviation from 100% is due to the
 # combination of self-shading and the fact that being at a tilt blocks off
-# the portion of the sky behind the row. The first effect is modeled with
+# the portion of the sky behind the row. Following the approach detailed in
+# [1]_, the first effect would be modeled with
 # :py:func:`pvlib.shading.sky_diffuse_passias` and the second with
 # :py:func:`pvlib.irradiance.isotropic`.
 
 plt.figure()
 for k in [1, 1.5, 2, 10]:
-    gcr = 1/k
+    gcr = 1 / k
     psi = shading.masking_angle_passias(surface_tilt, gcr)
     shading_loss = shading.sky_diffuse_passias(psi)
     transposition_ratio = irradiance.isotropic(surface_tilt, dhi=1.0)
-    relative_diffuse = transposition_ratio * (1-shading_loss) * 100  # %
+    relative_diffuse = transposition_ratio * (1 - shading_loss) * 100  # %
     plt.plot(surface_tilt, relative_diffuse, label=f'k={k}')
 
 plt.xlabel('Inclination angle [degrees]')
@@ -82,3 +89,37 @@
 # %%
 # As ``k`` decreases, GCR increases, so self-shading loss increases and
 # collected diffuse irradiance decreases.
+#
+# However, the pvlib-python authors believe that this approach is incorrect.
+#
+# Instead, the combination of inter-row shading from the previous row and the
+# surface tilt blocking the portion of the sky behind the row is obtained by
+# applying :py:func:`pvlib.shading.sky_diffuse_passias` on the sum of the
+# masking and surface tilt angles (see dashed curves in below figure). The
+# difference with the above approach is marginal for a ground coverage ratio
+# of 10%, but becomes very significant for high ground coverage ratios.
+#
+# Alternatively, one can also use :py:func:`bifacial.utils.vf_row_sky_2d_integ`
+# (see dotted curve in below figure), with very similar results except for the
+# highest ground coverage ratio. It is believed that the deviation is a result
+# of an approximation in :py:func:`pvlib.shading.masking_angle_passias` and
+# that :py:func:`bifacial.utils.vf_row_sky_2d_integ` provides the most accurate
+# result.
+
+color_cycler = cycler('color', ['blue', 'orange', 'green', 'red'])
+linestyle_cycler = cycler('linestyle', ['--', ':'])
+plt.rc('axes', prop_cycle=color_cycler * linestyle_cycler)
+plt.figure()
+for k in [1, 1.5, 2, 10]:
+    gcr = 1 / k
+    psi = shading.masking_angle_passias(surface_tilt, gcr)
+    vf1 = (1 - shading.sky_diffuse_passias(surface_tilt + psi)) * 100  # %
+    vf2 = bifacial.utils.vf_row_sky_2d_integ(surface_tilt, gcr) * 100  # %
+    plt.plot(surface_tilt, vf1, label=f'k={k} passias corrected')
+    plt.plot(surface_tilt, vf2, label=f'k={k} vf_row_sky_2d_integ')
+
+plt.xlabel('Inclination angle [degrees]')
+plt.ylabel('Relative diffuse irradiance [%]')
+plt.ylim(0, 105)
+plt.legend()
+plt.show()
diff --git a/docs/sphinx/source/whatsnew/v0.13.2.rst b/docs/sphinx/source/whatsnew/v0.13.2.rst
index e12f7277e5..f51639b5b1 100644
--- a/docs/sphinx/source/whatsnew/v0.13.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.13.2.rst
@@ -10,6 +10,8 @@ Breaking Changes
 
 Deprecations
 ~~~~~~~~~~~~
+* Deprecate :py:func:`~pvlib.shading.sky_diffuse_passias` in favor of
+:py:func:`~pvlib.bifacial.utils.vf_row_sky_2d_integ`. (:pull:`2584`)
 
 
 Bug fixes
@@ -33,7 +35,7 @@ Enhancements
 Documentation
 ~~~~~~~~~~~~~
 * Provide an overview of single-diode modeling functionality in :ref:`singlediode`. (:pull:`2565`)
-
+* Update gallery example of modeling diffuse shading loss with correct usage. (:pull:`2584`)
 
 Testing
 ~~~~~~~
diff --git a/pvlib/shading.py b/pvlib/shading.py
index 42aef892f7..e7c56ff0e2 100644
--- a/pvlib/shading.py
+++ b/pvlib/shading.py
@@ -5,6 +5,8 @@
 
 import numpy as np
 import pandas as pd
+
+from pvlib._deprecation import deprecated
 from pvlib.tools import sind, cosd
 
 
@@ -126,6 +128,7 @@ def masking_angle_passias(surface_tilt, gcr):
     --------
     masking_angle
     sky_diffuse_passias
+    bifacial.utils.vf_row_sky_2d_integ
 
     Notes
     -----
@@ -191,6 +194,12 @@ def masking_angle_passias(surface_tilt, gcr):
     return np.degrees(psi_avg)
 
 
+@deprecated(
+    "0.13.2",
+    pending=True,
+    alternative="bifacial.utils.vf_row_sky_2d_integ",
+    addendum="See also the 'Diffuse Self-Shading' Example in the Gallery.",
+)
 def sky_diffuse_passias(masking_angle):
     r"""
     The diffuse irradiance loss caused by row-to-row sky diffuse shading.
@@ -198,10 +207,11 @@ def sky_diffuse_passias(masking_angle):
     Even when the sun is high in the sky, a row's view of the sky dome will
     be partially blocked by the row in front. This causes a reduction in the
     diffuse irradiance incident on the module. The reduction depends on the
-    masking angle, the elevation angle from a point on the shaded module to
-    the top of the shading row. In [1]_ the masking angle is calculated as
-    the average across the module height. SAM assumes the "worst-case" loss
-    where the masking angle is calculated for the bottom of the array [2]_.
+    masking angle, the surface tilt and the elevation angle from a point on
+    the shaded module to the top of the shading row. In [1]_ the masking
+    angle is calculated as the average across the module height. SAM assumes
+    the "worst-case" loss where the masking angle is calculated for the
+    bottom of the array [2]_.
 
     This function, as in [1]_, makes the assumption that sky diffuse
     irradiance is isotropic.
@@ -209,8 +219,8 @@ def sky_diffuse_passias(masking_angle):
     Parameters
     ----------
     masking_angle : numeric
-        The elevation angle below which diffuse irradiance is blocked
-        [degrees].
+        The sum of the surface tilt and the elevation angle below which
+        diffuse irradiance is blocked [degrees].
 
     Returns
     -------
@@ -221,6 +231,7 @@ def sky_diffuse_passias(masking_angle):
     --------
     masking_angle
     masking_angle_passias
+    bifacial.utils.vf_row_sky_2d_integ
 
     References
     ----------