diff --git a/test/grid/grid/test_areas.py b/test/grid/grid/test_areas.py
index bce897a25..f11cbde8d 100644
--- a/test/grid/grid/test_areas.py
+++ b/test/grid/grid/test_areas.py
@@ -116,7 +116,6 @@ def test_latlon_bounds_populate_bounds_MPAS(gridpath):
     uxgrid = ux.open_grid(gridpath("mpas", "QU", "oQU480.231010.nc"))
     bounds_xarray = uxgrid.bounds
 
-
 def _sum_quadrature_jacobians(x, y, z, quadrature_rule, order):
     """Independently sum the Jacobian at each quadrature point of a triangle.
 
@@ -163,3 +162,10 @@ def test_calculate_face_area_jacobian_is_quadrature_sum(quadrature_rule, order):
 
     assert jacobian > 0
     nt.assert_allclose(jacobian, expected, rtol=1e-12)
+
+
+def test_face_areas_fesom(gridpath):
+    """Ensure correct total area for FESOM grid (~8.3780 sr). Regression test for #425."""
+    uxgrid = ux.open_grid(gridpath("ugrid", "fesom", "fesom.mesh.diag.nc"))
+    total_area = uxgrid.calculate_total_face_area()
+    nt.assert_almost_equal(total_area, 8.3780, decimal=4)
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 6fa02b069..9bd17e8a2 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -1958,7 +1958,12 @@ def calculate_total_face_area(
         """Calculate the total surface area of all the faces in a mesh.
 
         Equivalent to ``self.compute_face_areas(...).sum()``; provided as a
-        convenience.
+        convenience. (Note: for HEALPix grids, when called with default arguments,
+        this method actually returns ``self.face_areas.sum()`` instead,
+        which respects HEALPix equal-area property.)
+
+        Additionally, raises a warning if the result is larger than
+        the total area of a sphere (4 * pi * self.sphere_radius**2).
 
         Parameters
         ----------
@@ -1986,9 +1991,9 @@ def calculate_total_face_area(
             and order == 4
             and not latitude_adjusted_area
         ):
-            return np.sum(self.face_areas.values)
+            result = np.sum(self.face_areas.values)
 
-        return np.sum(
+        result = np.sum(
             self.compute_face_areas(
                 quadrature_rule=quadrature_rule,
                 order=order,
@@ -1996,6 +2001,15 @@ def calculate_total_face_area(
             )
         )
 
+        RTOL = 1e-6  # 1e-7 had warnings in existing CI tests (as of 2026-08-05), 1e-6 did not.
+        if result > 4 * np.pi * self.sphere_radius**2 * (1 + RTOL):
+            warnings.warn(
+                f"Total face area (={result}) exceeds the surface area of the whole sphere "
+                f"(={4 * np.pi * self.sphere_radius**2}) (with sphere_radius={self.sphere_radius}).",
+            )
+
+        return result
+
     def compute_face_areas(
         self,
         quadrature_rule: str = "triangular",