diff --git a/src/meshes/CutCellMeshes/layered.jl b/src/meshes/CutCellMeshes/layered.jl
index 7e5a96897..ac277c1f1 100644
--- a/src/meshes/CutCellMeshes/layered.jl
+++ b/src/meshes/CutCellMeshes/layered.jl
@@ -89,7 +89,7 @@ function layered_mesh(
     mesh::UnstructuredMesh{3},
     surfaces::Vector{<:PolygonalSurface};
     min_cut_fraction::Real = 0.01,
-    merge_cells::Bool = true
+    merge_cells::Bool = false
 )
     T = Float64
     n_surfaces = length(surfaces)
diff --git a/src/meshes/CutCellMeshes/merge_faces.jl b/src/meshes/CutCellMeshes/merge_faces.jl
index 22a21e6e9..532f4186e 100644
--- a/src/meshes/CutCellMeshes/merge_faces.jl
+++ b/src/meshes/CutCellMeshes/merge_faces.jl
@@ -203,11 +203,24 @@ function _merge_face_group!(
             end
         end
 
-        # Extract the boundary polygon by chaining boundary edges
+        # Check that the boundary forms a simple polygon: every boundary
+        # node must have exactly two boundary edges.  If any node has a
+        # different degree the merged outline is non-simple and cannot be
+        # properly chained into a single polygon.
+        bnd_deg = Dict{Int, Int}()
+        for (a, b) in boundary_edges
+            bnd_deg[a] = get(bnd_deg, a, 0) + 1
+            bnd_deg[b] = get(bnd_deg, b, 0) + 1
+        end
+        if any(d != 2 for d in values(bnd_deg))
+            continue  # non-simple boundary; leave faces un-merged
+        end
+
+        # Extract the boundary polygon by chaining boundary edges.
         merged_nodes = _chain_boundary_edges(boundary_edges)
-        if isempty(merged_nodes)
-            # Fallback: use all nodes sorted by angle
-            merged_nodes = collect(all_nodes)
+        if isempty(merged_nodes) || length(merged_nodes) != length(bnd_deg)
+            # Chain did not cover all boundary nodes
+            continue  # leave faces un-merged
         end
 
         # Order the polygon nodes correctly
diff --git a/test/cut_cell_meshes.jl b/test/cut_cell_meshes.jl
index a04f44277..837fdc6f3 100644
--- a/test/cut_cell_meshes.jl
+++ b/test/cut_cell_meshes.jl
@@ -848,6 +848,43 @@ import Jutul.CutCellMeshes: PlaneCut, PolygonalSurface, cut_mesh, layered_mesh,
         end
     end
 
+    @testset "merge_coplanar_faces after layered_mesh preserves volume" begin
+        import Jutul.CutCellMeshes: merge_coplanar_faces
+
+        # Reproduce a bug where merge_coplanar_faces created holes
+        # when applied after layered_mesh with a non-planar surface.  The root
+        # cause was non-simple polygon boundaries (nodes with degree != 2) that
+        # the greedy chain algorithm could not trace correctly.
+        L = 1000.0
+        g = CartesianMesh((3, 3, 3), (L, L, 100.0))
+        mesh = UnstructuredMesh(g)
+
+        N = 4
+        xs = range(0.0, L, length = N)
+        ys = range(0.0, L, length = N)
+
+        # Use a deterministic non-planar surface that exercises the fix
+        zvals = [50.0 + 5.0 * sin(2π * (i - 1) / (N - 1)) * cos(2π * (j - 1) / (N - 1))
+                 for i in 1:N, j in 1:N]
+        surf = depth_grid_to_surface(xs, ys, zvals)
+        result, _ = layered_mesh(mesh, [surf])
+
+        geo_before = tpfv_geometry(result)
+        vol_before = sum(geo_before.volumes)
+
+        merged = merge_coplanar_faces(result)
+        geo_after = tpfv_geometry(merged)
+        vol_after = sum(geo_after.volumes)
+
+        # Volume must be preserved (tight tolerance)
+        @test vol_after ≈ vol_before rtol = 1e-6
+        # All cell volumes must remain positive
+        @test all(geo_after.volumes .> 0)
+        # Face count must not increase
+        @test number_of_faces(merged) <= number_of_faces(result)
+        @test number_of_boundary_faces(merged) <= number_of_boundary_faces(result)
+    end
+
     @testset "merge_faces with extra_out" begin
         g = CartesianMesh((3, 3, 3))
         mesh = UnstructuredMesh(g)