diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst
index fff7ba5d46..345ffd53c1 100644
--- a/docs/sphinx/source/whatsnew/v0.16.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.2.rst
@@ -16,6 +16,8 @@ Bug fixes
 ~~~~~~~~~
 * :py:func:`pvlib.tracking.calc_surface_orientation` returns scalar when input
   is scalar (:issue:`2881`, :pull:`2882`)
+* Fix :py:func:`pvlib.solarposition.solar_azimuth_analytical` returning south
+  instead of north at solar noon.
 
 Enhancements
 ~~~~~~~~~~~~
@@ -62,6 +64,7 @@ Contributors
 * Chen Yi (:ghuser:`chenyi6758-source`)
 * Rajiv Daxini (:ghuser:`rdaxini`)
 * jason-writer (:ghuser:`jason-writer`)
+* :ghuser:`vera-rubin`
 
 
 
diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py
index 4441ed01f1..6634c8c7a2 100644
--- a/pvlib/solarposition.py
+++ b/pvlib/solarposition.py
@@ -1308,6 +1308,9 @@ def solar_azimuth_analytical(latitude, hourangle, declination, zenith):
     with np.errstate(invalid='ignore'):
         sign_ha = np.sign(hourangle)
 
+    # At solar noon, use the morning branch so that north is zero radians.
+    sign_ha = sign_ha - 1.0 * (sign_ha == 0)
+
     return sign_ha * np.arccos(cos_azi) + np.pi
 
 
diff --git a/tests/test_solarposition.py b/tests/test_solarposition.py
index 9dc0488fe5..824e046e28 100644
--- a/tests/test_solarposition.py
+++ b/tests/test_solarposition.py
@@ -718,6 +718,22 @@ def test_analytical_azimuth():
     assert not np.isnan(azimuths).any()
 
 
+@pytest.mark.parametrize('latitude,declination,expected', [
+    (-35., 0., 0.),
+    (35., 0., 180.),
+    (10., 20., 0.),
+])
+@pytest.mark.parametrize('hourangle', [
+    0., -0., np.array([0.]), np.array([-0.]), pd.Series([0., -0.]),
+])
+def test_analytical_azimuth_noon(latitude, declination, expected, hourangle):
+    # At noon, the sun is due north if declination exceeds latitude.
+    zenith = np.radians(abs(latitude - declination))
+    azimuth = solarposition.solar_azimuth_analytical(
+        np.radians(latitude), hourangle, np.radians(declination), zenith)
+    assert_allclose(azimuth, np.radians(expected), atol=1e-8)
+
+
 def test_hour_angle():
     """
     Test conversion from hours to hour angles in degrees given the following