diff --git a/pyproject.toml b/pyproject.toml
index 0a2c87a..be18ea2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -44,6 +44,7 @@ dependencies = [
     "ipython >= 8.0.0, < 10",
     "lib3mf >= 2.4.1",
     "ocpsvg >= 0.5, < 0.6",
+    "ocp_gordon >= 0.1.10",
     "trianglesolver",
     "sympy",
     "scipy",
diff --git a/src/build123d/topology/two_d.py b/src/build123d/topology/two_d.py
index 306c5b8..d265114 100644
--- a/src/build123d/topology/two_d.py
+++ b/src/build123d/topology/two_d.py
@@ -105,6 +105,7 @@ from OCP.TopExp import TopExp
 from OCP.TopoDS import TopoDS, TopoDS_Face, TopoDS_Shape, TopoDS_Shell, TopoDS_Solid
 from OCP.TopTools import TopTools_IndexedDataMapOfShapeListOfShape, TopTools_ListOfShape
 from typing_extensions import Self
+from ocp_gordon import interpolate_curve_network
 
 from build123d.build_enums import (
     CenterOf,
@@ -864,6 +865,38 @@ class Face(Mixin2D, Shape[TopoDS_Face]):
             raise ValueError("Can't extrude empty object")
         return Face(TopoDS.Face_s(_extrude_topods_shape(obj.wrapped, direction)))
 
+    @classmethod
+    def gordon_surface(
+        cls,
+        profiles: Iterable[Edge],
+        guides: Iterable[Edge],
+        tolerance: float = 3e-4,
+    ) -> Face:
+        """
+        Creates a Gordon surface from a network of profile and guide curves.
+
+        Args:
+            profiles (Iterable[Edge]): Edges representing profile curves.
+            guides (Iterable[Edge]): Edges representing guide curves.
+            tolerance (float, optional): Tolerance for surface creation and
+                intersection calculations.
+
+        Returns:
+            Face: the interpolated Gordon surface
+        """
+        ocp_profiles = [BRep_Tool.Curve_s(edge.wrapped, 0, 1) for edge in profiles]
+        ocp_guides = [BRep_Tool.Curve_s(edge.wrapped, 0, 1) for edge in guides]
+
+        gordon_bspline_surface = interpolate_curve_network(
+            ocp_profiles, ocp_guides, tolerance=tolerance
+        )
+
+        return cls(
+            BRepBuilderAPI_MakeFace(
+                gordon_bspline_surface, Precision.Confusion_s()
+            ).Face()
+        )
+
     @classmethod
     def make_bezier_surface(
         cls,
diff --git a/tests/test_direct_api/test_face.py b/tests/test_direct_api/test_face.py
index 82460c0..3e14740 100644
--- a/tests/test_direct_api/test_face.py
+++ b/tests/test_direct_api/test_face.py
@@ -502,6 +502,78 @@ class TestFace(unittest.TestCase):
                 ]
             )
 
+    def test_gordon_surface(self):
+        def create_test_curves(
+            num_profiles: int = 3,
+            num_guides: int = 4,
+            u_range: float = 1.0,
+            v_range: float = 1.0,
+        ):
+            profiles: list[Edge] = []
+            guides: list[Edge] = []
+
+            intersection_points = [
+                [(0.0, 0.0, 0.0) for _ in range(num_guides)]
+                for _ in range(num_profiles)
+            ]
+
+            for i in range(num_profiles):
+                for j in range(num_guides):
+                    u = i * u_range / (num_profiles - 1)
+                    v = j * v_range / (num_guides - 1)
+                    z = 0.2 * math.sin(u * math.pi) * math.cos(v * math.pi)
+                    intersection_points[i][j] = (u, v, z)
+
+            for i in range(num_profiles):
+                points = [intersection_points[i][j] for j in range(num_guides)]
+                profiles.append(Spline(points))
+
+            for j in range(num_guides):
+                points = [intersection_points[i][j] for i in range(num_profiles)]
+                guides.append(Spline(points))
+
+            return profiles, guides
+
+        profiles, guides = create_test_curves()
+
+        tolerance = 3e-4
+        gordon_surface = Face.gordon_surface(profiles, guides, tolerance=tolerance)
+
+        self.assertIsInstance(
+            gordon_surface, Face, "The returned object should be a Face."
+        )
+
+        def point_at_uv_against_expected(u: float, v: float, expected_point: Vector):
+            point_at_uv = gordon_surface.position_at(u, v)
+            self.assertAlmostEqual(
+                point_at_uv.X,
+                expected_point.X,
+                delta=tolerance,
+                msg=f"X coordinate mismatch at ({u},{v})",
+            )
+            self.assertAlmostEqual(
+                point_at_uv.Y,
+                expected_point.Y,
+                delta=tolerance,
+                msg=f"Y coordinate mismatch at ({u},{v})",
+            )
+            self.assertAlmostEqual(
+                point_at_uv.Z,
+                expected_point.Z,
+                delta=tolerance,
+                msg=f"Z coordinate mismatch at ({u},{v})",
+            )
+
+        point_at_uv_against_expected(
+            u=1.0, v=0.0, expected_point=profiles[0].position_at(1.0)
+        )
+        point_at_uv_against_expected(
+            u=0.0, v=1.0, expected_point=guides[0].position_at(1.0)
+        )
+        point_at_uv_against_expected(
+            u=1.0, v=1.0, expected_point=profiles[-1].position_at(1.0)
+        )
+
     # def test_to_arcs(self):
     #     with BuildSketch() as bs:
     #         with BuildLine() as bl: