diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst
index 64132a2547..a2eebbac2a 100644
--- a/docs/sphinx/source/whatsnew/v0.16.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.2.rst
@@ -21,6 +21,11 @@ Bug fixes
   is scalar (:issue:`2881`, :pull:`2882`)
 * :py:func:`pvlib.scaling.wvm` no longer returns all NaNs when ``positions``
   contains a single site (:issue:`2895`, :pull:`2896`)
+* :py:func:`pvlib.irradiance.beam_component` and
+  :py:func:`pvlib.irradiance.poa_components` now clip the angle of incidence
+  projection before multiplying by ``dni``, so negative DNI no longer yields
+  positive beam irradiance when the sun is behind the surface.
+  (:issue:`2903`, :pull:`2913`)
 
 Enhancements
 ~~~~~~~~~~~~
@@ -69,3 +74,4 @@ Contributors
 * jason-writer (:ghuser:`jason-writer`)
 * :ghuser:`bfs2021`
 * :ghuser:`xiajiahang`
+* :ghuser:`Noumenon-ai`
diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index 38bf122b91..d77c7b5d1d 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -268,9 +268,9 @@ def beam_component(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
     beam : numeric
         Beam component. [Wm⁻²]
     """
-    beam = dni * aoi_projection(surface_tilt, surface_azimuth,
+    projection = aoi_projection(surface_tilt, surface_azimuth,
                                 solar_zenith, solar_azimuth)
-    beam = np.maximum(beam, 0)
+    beam = dni * np.maximum(projection, 0)
 
     return beam
 
@@ -545,7 +545,7 @@ def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
 
     Notes
     ------
-    Negative beam irradiation due to AOI > 90° or AOI < 0° is set to zero.
+    Beam irradiance is set to zero when the sun is behind the surface.
     '''
 
     if isinstance(poa_sky_diffuse, dict):
@@ -558,7 +558,7 @@ def poa_components(aoi, dni, poa_sky_diffuse, poa_ground_diffuse):
         sky_components = {}
         total_poa_sky_diffuse = poa_sky_diffuse
 
-    poa_direct = np.maximum(dni * np.cos(np.radians(aoi)), 0)
+    poa_direct = dni * np.maximum(np.cos(np.radians(aoi)), 0)
     poa_diffuse = total_poa_sky_diffuse + poa_ground_diffuse
     poa_global = poa_direct + poa_diffuse
 
diff --git a/tests/test_irradiance.py b/tests/test_irradiance.py
index b0d9863ad5..b2d5f4ed92 100644
--- a/tests/test_irradiance.py
+++ b/tests/test_irradiance.py
@@ -693,6 +693,25 @@ def test_get_total_irradiance_missing_airmass():
                                   'poa_ground_diffuse']
 
 
+@pytest.mark.parametrize('array_type', [np.array, pd.Series])
+def test_beam_component_negative_dni(array_type):
+    solar_zenith = array_type([120., 120., 60., 60., 90., 90., np.nan])
+    dni = array_type([-10., 100., -10., 100., -10., 100., 100.])
+    result = irradiance.beam_component(0, 180, solar_zenith, 180, dni)
+    assert_allclose(result, [0., 0., -5., 50., 0., 0., np.nan], atol=1e-12)
+
+
+@pytest.mark.parametrize('array_type', [np.array, pd.Series])
+def test_poa_components_negative_dni(array_type):
+    aoi = array_type([120., 120., 60., 60., 90., 90., np.nan])
+    dni = array_type([-10., 100., -10., 100., -10., 100., 100.])
+    result = irradiance.poa_components(aoi, dni, 10., 5.)
+    assert_allclose(result['poa_direct'],
+                    [0., 0., -5., 50., 0., 0., np.nan], atol=1e-12)
+    assert_allclose(result['poa_global'],
+                    [15., 15., 10., 65., 15., 15., np.nan])
+
+
 def test_poa_components(irrad_data, ephem_data, dni_et, relative_airmass):
     aoi = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
                          ephem_data['azimuth'])