diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py
index 9931889c..bc74dac1 100644
--- a/boule/tests/test_ellipsoid.py
+++ b/boule/tests/test_ellipsoid.py
@@ -8,6 +8,7 @@
 import numpy.testing as npt
 
 from .. import Ellipsoid, ELLIPSOIDS
+from .utils import normal_gravity_surface
 
 
 ELLIPSOID_NAMES = [e.name for e in ELLIPSOIDS]
@@ -239,3 +240,15 @@ def test_geocentric_radius_geocentric_pole_equator(ellipsoid):
     npt.assert_allclose(
         radius_true, ellipsoid.geocentric_radius(latitude, geodetic=False)
     )
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_normal_gravity_against_somigliana(ellipsoid):
+    """
+    Check if normal gravity on the surface satisfies Somigliana equation
+    """
+    latitude = np.linspace(-90, 90, 181)
+    npt.assert_allclose(
+        ellipsoid.normal_gravity(latitude, height=0),
+        normal_gravity_surface(latitude, ellipsoid),
+    )
diff --git a/boule/tests/utils.py b/boule/tests/utils.py
new file mode 100644
index 00000000..335aff5d
--- /dev/null
+++ b/boule/tests/utils.py
@@ -0,0 +1,24 @@
+"""
+Shared utility functions for testing.
+"""
+import numpy as np
+
+
+def normal_gravity_surface(latitude, ellipsoid):
+    """
+    Computes normal gravity on the surface of the ellipsoid [mGal]
+
+    Uses the closed-form Somigliana equation [Hofmann-WellenhofMoritz2006]_.
+    """
+    latitude_radians = np.radians(latitude)
+    coslat = np.cos(latitude_radians)
+    sinlat = np.sin(latitude_radians)
+    gravity = (
+        ellipsoid.semimajor_axis * ellipsoid.gravity_equator * coslat ** 2
+        + ellipsoid.semiminor_axis * ellipsoid.gravity_pole * sinlat ** 2
+    ) / np.sqrt(
+        ellipsoid.semimajor_axis ** 2 * coslat ** 2
+        + ellipsoid.semiminor_axis ** 2 * sinlat ** 2
+    )
+    # Convert to mGal
+    return 1e5 * gravity