diff --git a/docs/sphinx/source/whatsnew/v0.16.0.rst b/docs/sphinx/source/whatsnew/v0.16.0.rst
index 6d62b9e1c6..9bbafb16ca 100644
--- a/docs/sphinx/source/whatsnew/v0.16.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.0.rst
@@ -19,6 +19,9 @@ Deprecations
 
 Bug fixes
 ~~~~~~~~~
+* :py:func:`~pvlib.irradiance.perez` no longer raises ``ZeroDivisionError`` on
+  scalar ``dhi=0`` input (e.g. nighttime); the scalar path now returns the same
+  finite value as the array path. (:pull:`2826`)
 
 
 Enhancements
@@ -48,3 +51,4 @@ Maintenance
 Contributors
 ~~~~~~~~~~~~
 * Carolina Crespo (:ghuser:`cbcrespo`)
+* Andrew Chen (:ghuser:`chuenchen309`)
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 98bd948286..e10da59609 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -1146,8 +1146,12 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     delta = dhi * airmass / dni_extra
 
     # epsilon is the sky's "clearness"
-    with np.errstate(invalid='ignore'):
-        eps = ((dhi + dni) / dhi + kappa * (z ** 3)) / (1 + kappa * (z ** 3))
+    # np.true_divide so a Python-scalar dhi=0 yields inf/nan like the array
+    # path (handled below via digitize) instead of raising ZeroDivisionError,
+    # which np.errstate cannot suppress for native scalar division.
+    with np.errstate(invalid='ignore', divide='ignore'):
+        eps = (np.true_divide(dhi + dni, dhi) + kappa * (z ** 3)) \
+            / (1 + kappa * (z ** 3))
 
     # numpy indexing below will not work with a Series
     if isinstance(eps, pd.Series):
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index fbe9906358..0e33f60f73 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -429,6 +429,16 @@ def test_perez_scalar():
     assert_allclose(out, 109.084332)
 
 
+def test_perez_scalar_dhi_zero():
+    # dhi=0 (e.g. nighttime) is documented-valid input (dhi >= 0). A native
+    # Python scalar division raised ZeroDivisionError (which np.errstate cannot
+    # suppress) while the array path returns a finite value; they must match.
+    out = irradiance.perez(30, 180, 0.0, 800.0, 1400.0, 40.0, 120.0, 1.5)
+    expected = irradiance.perez(30, 180, np.array([0.0]), np.array([800.0]),
+                                1400.0, 40.0, 120.0, np.array([1.5]))
+    assert_allclose(out, expected[0])
+
+
 def test_perez_driesse_scalar():
     # copied values from fixtures
     out = irradiance.perez_driesse(40, 180, 118.458, 939.954,