diff --git a/geoapps_utils/utils/transformations.py b/geoapps_utils/utils/transformations.py
index c85429af..af3146fc 100644
--- a/geoapps_utils/utils/transformations.py
+++ b/geoapps_utils/utils/transformations.py
@@ -10,6 +10,7 @@
 
 import numpy as np
 import scipy.sparse as ssp
+from geoh5py.objects import Surface
 
 
 def z_rotation_matrix(angle: np.ndarray) -> ssp.csr_matrix:
@@ -195,3 +196,24 @@ def cartesian_normal_to_direction_and_dip(normals: np.ndarray) -> np.ndarray:
     direction_and_dip = spherical_normal_to_direction_and_dip(spherical_normals[:, 1:])
 
     return direction_and_dip
+
+
+def compute_normals(surface: Surface) -> np.ndarray:
+    """
+    Compute normals for each triangle in a surface.
+
+    :param surface: Surface object containing vertices and cells.
+
+    :returns: Array of shape (n, 3) representing the normals for
+        each cell in the mesh.
+    """
+
+    vertices = surface.vertices
+    cells = surface.cells
+
+    v1 = vertices[cells[:, 1]] - vertices[cells[:, 0]]
+    v2 = vertices[cells[:, 2]] - vertices[cells[:, 0]]
+    normals = np.cross(v1, v2)
+    normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis]
+
+    return normals
diff --git a/tests/transformations_test.py b/tests/transformations_test.py
index 2f4c7ec0..7a4b62f4 100644
--- a/tests/transformations_test.py
+++ b/tests/transformations_test.py
@@ -14,11 +14,12 @@
 import pytest
 from geoh5py import Workspace
 from geoh5py.groups.property_group import GroupTypeEnum
-from geoh5py.objects import Points
+from geoh5py.objects import Points, Surface
 
 from geoapps_utils.utils.transformations import (
     cartesian_normal_to_direction_and_dip,
     cartesian_to_spherical,
+    compute_normals,
     rotate_xyz,
     spherical_normal_to_direction_and_dip,
     x_rotation_matrix,
@@ -143,3 +144,37 @@ def test_spherical_values(tmp_path):  # pylint: disable=too-many-locals
             properties=[direction, dip],
             property_group_type=GroupTypeEnum.DIPDIR,
         )
+
+
+def create_surface(workspace, upside_down=False):
+    k = -5.0 if upside_down else 5.0
+    vertices = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [10.0, 0.0, 0.0],
+            [10.0, 10.0, 0.0],
+            [0.0, 10.0, 0.0],
+            [5.0, 5.0, k],
+        ]
+    )
+    triangles = np.array(
+        [
+            [0, 1, 4],
+            [1, 2, 4],
+            [2, 3, 4],
+            [3, 0, 4],
+        ]
+    )
+    surface = Surface.create(workspace, vertices=vertices, cells=triangles)
+
+    return surface
+
+
+def test_compute_normals(tmp_path):
+    ws = Workspace(tmp_path / "test.geoh5")
+    surface = create_surface(ws, upside_down=False)
+    normals = compute_normals(surface)
+
+    h = np.sqrt(2) / 2
+    validation = np.array([[0, -h, h], [h, 0, h], [0, h, h], [-h, 0, h]])
+    assert np.allclose(normals, validation)