diff --git a/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py b/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py
new file mode 100644
index 0000000000..5525c49453
--- /dev/null
+++ b/docs/examples/bifacial/plot_pvfactors_fixed_tilt.py
@@ -0,0 +1,68 @@
+"""
+Fixed-Tilt Simulation with pvfactors
+====================================
+
+Modeling the irradiance on the rear side of a fixed-tilt array.
+"""
+
+# %%
+# Because pvfactors was originally designed for modeling single-axis
+# tracking systems, it's not necessarily obvious how to use it to model
+# fixed-tilt systems correctly.
+# This example shows how to model rear-side irradiance on a fixed-tilt
+# array using :py:func:`pvlib.bifacial.pvfactors.pvfactors_timeseries`.
+
+import pandas as pd
+from pvlib import location
+from pvlib.bifacial.pvfactors import pvfactors_timeseries
+import matplotlib.pyplot as plt
+import warnings
+
+# supressing shapely warnings that occur on import of pvfactors
+warnings.filterwarnings(action='ignore', module='pvfactors')
+
+# %%
+# First, generate the usual modeling inputs:
+
+times = pd.date_range('2021-06-21', '2021-06-22', freq='1T', tz='Etc/GMT+5')
+loc = location.Location(latitude=40, longitude=-80, tz=times.tz)
+sp = loc.get_solarposition(times)
+cs = loc.get_clearsky(times)
+
+# example array geometry
+pvrow_height = 1
+pvrow_width = 4
+pitch = 10
+gcr = pvrow_width / pitch
+axis_azimuth = 180
+albedo = 0.2
+
+# %%
+# Now the trick: since pvfactors only wants to model single-axis tracking
+# arrays, we have to pretend our fixed tilt array is a single-axis tracking
+# array that never rotates.  In that case, the "axis of rotation" is
+# along the length of the row, with ``axis_azimuth`` 90 degrees offset from the
+# fixed ``surface_azimuth``.
+
+irrad = pvfactors_timeseries(
+    solar_azimuth=sp['azimuth'],
+    solar_zenith=sp['apparent_zenith'],
+    surface_azimuth=180,  # south-facing array
+    surface_tilt=20,
+    axis_azimuth=90,  # 90 degrees off from surface_azimuth.  270 is ok too
+    timestamps=times,
+    dni=cs['dni'],
+    dhi=cs['dhi'],
+    gcr=gcr,
+    pvrow_height=pvrow_height,
+    pvrow_width=pvrow_width,
+    albedo=albedo,
+    n_pvrows=3,
+    index_observed_pvrow=1
+)
+
+# turn into pandas DataFrame
+irrad = pd.concat(irrad, axis=1)
+
+irrad[['total_inc_back', 'total_abs_back']].plot()
+plt.ylabel('Irradiance [W m$^{-2}$]')
diff --git a/docs/sphinx/source/whatsnew/v0.9.2.rst b/docs/sphinx/source/whatsnew/v0.9.2.rst
index 893cc65639..44ca278fb2 100644
--- a/docs/sphinx/source/whatsnew/v0.9.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.9.2.rst
@@ -27,6 +27,9 @@ Testing
 
 Documentation
 ~~~~~~~~~~~~~
+* Add gallery example of simulating rearside irradiance for a fixed-tilt
+  array with pvfactors (:pull:`1470`)
+
 
 Benchmarking
 ~~~~~~~~~~~~~
diff --git a/pvlib/bifacial/pvfactors.py b/pvlib/bifacial/pvfactors.py
index 3c70b779fb..ed09879a21 100644
--- a/pvlib/bifacial/pvfactors.py
+++ b/pvlib/bifacial/pvfactors.py
@@ -35,6 +35,8 @@ def pvfactors_timeseries(
     axis_azimuth: float
         Azimuth angle of the rotation axis of the PV modules, using pvlib's
         convention (deg). This is supposed to be fixed for all timestamps.
+        When modeling fixed-tilt arrays, set this value to be 90 degrees
+        clockwise from ``surface_azimuth``.
     timestamps: datetime or DatetimeIndex
         List of simulation timestamps
     dni: numeric