diff --git a/src/NamedGraphs.jl b/src/NamedGraphs.jl
index e5985fb..8de2717 100644
--- a/src/NamedGraphs.jl
+++ b/src/NamedGraphs.jl
@@ -113,6 +113,7 @@ include("distance.jl")
 include("distances_and_capacities.jl")
 include(joinpath("steiner_tree", "steiner_tree.jl"))
 include(joinpath("traversals", "dfs.jl"))
+include(joinpath("traversals", "trees_and_forests.jl"))
 include("namedgraph.jl")
 include(joinpath("generators", "named_staticgraphs.jl"))
 include(joinpath("Graphs", "generators", "staticgraphs.jl"))
diff --git a/src/traversals/trees_and_forests.jl b/src/traversals/trees_and_forests.jl
new file mode 100644
index 0000000..de433ee
--- /dev/null
+++ b/src/traversals/trees_and_forests.jl
@@ -0,0 +1,56 @@
+abstract type SpanningTreeAlgorithm end
+
+struct BFS <: SpanningTreeAlgorithm end
+struct RandomBFS <: SpanningTreeAlgorithm end
+struct DFS <: SpanningTreeAlgorithm end
+
+default_spanning_tree_alg() = BFS()
+
+default_root_vertex(g) = last(findmax(eccentricities(g)))
+
+function spanning_tree(
+  g::AbstractNamedGraph; alg=default_spanning_tree_alg(), root_vertex=default_root_vertex(g)
+)
+  return spanning_tree(alg, g; root_vertex)
+end
+
+@traitfn function spanning_tree(
+  ::BFS, g::AbstractNamedGraph::(!IsDirected); root_vertex=default_root_vertex(g)
+)
+  return undirected_graph(bfs_tree(g, root_vertex))
+end
+
+@traitfn function spanning_tree(
+  ::RandomBFS, g::AbstractNamedGraph::(!IsDirected); root_vertex=default_root_vertex(g)
+)
+  return undirected_graph(random_bfs_tree(g, root_vertex))
+end
+
+@traitfn function spanning_tree(
+  ::DFS, g::AbstractNamedGraph::(!IsDirected); root_vertex=default_root_vertex(g)
+)
+  return undirected_graph(dfs_tree(g, root_vertex))
+end
+
+#Given a graph, split it into its connected components, construct a spanning tree, using the function spanning_tree, over each of them
+# and take the union.
+function spanning_forest(g::AbstractNamedGraph; spanning_tree=spanning_tree)
+  return reduce(union, (spanning_tree(g[vs]) for vs in connected_components(g)))
+end
+
+#Given an undirected graph g with vertex set V, build a set of forests (each with vertex set V) which covers all edges in g
+# (see https://en.wikipedia.org/wiki/Arboricity) We do not find the minimum but our tests show this algorithm performs well
+function forest_cover(g::AbstractNamedGraph; spanning_tree=spanning_tree)
+  edges_collected = edgetype(g)[]
+  remaining_edges = edges(g)
+  forests = NamedGraph[]
+  while !isempty(remaining_edges)
+    g_reduced = rem_edges(g, edges_collected)
+    g_reduced_spanning_forest = spanning_forest(g_reduced; spanning_tree)
+    push!(edges_collected, edges(g_reduced_spanning_forest)...)
+    push!(forests, g_reduced_spanning_forest)
+    setdiff!(remaining_edges, edges(g_reduced_spanning_forest))
+  end
+
+  return forests
+end
diff --git a/test/test_namedgraph.jl b/test/test_namedgraph.jl
index b52d9a8..7885f05 100644
--- a/test/test_namedgraph.jl
+++ b/test/test_namedgraph.jl
@@ -390,17 +390,16 @@ end
     g = NamedGraph(path_graph(4), ["A", "B", "C", "D"])
 
     part1, part2, flow = GraphsFlows.mincut(g, "A", "D")
-    @test issetequal(part1, ["A"])
-    @test issetequal(part2, ["B", "C", "D"])
+    @test "A" ∈ part1
+    @test "D" ∈ part2
     @test flow == 1
 
     part1, part2 = mincut_partitions(g, "A", "D")
-    @test issetequal(part1, ["A"])
-    @test issetequal(part2, ["B", "C", "D"])
+    @test "A" ∈ part1
+    @test "D" ∈ part2
 
     part1, part2 = mincut_partitions(g)
-    @test issetequal(part1, ["B", "C", "D"])
-    @test issetequal(part2, ["A"])
+    @test issetequal(vcat(part1, part2), vertices(g))
 
     weights_dict = Dict{Tuple{String,String},Float64}()
     weights_dict["A", "B"] = 3
@@ -411,17 +410,17 @@ end
 
     for weights in (weights_dict, weights_dictionary)
       part1, part2, flow = GraphsFlows.mincut(g, "A", "D", weights)
-      @test issetequal(part1, ["A", "B"])
-      @test issetequal(part2, ["C", "D"])
+      @test issetequal(part1, ["A", "B"]) || issetequal(part1, ["C", "D"])
+      @test issetequal(vcat(part1, part2), vertices(g))
       @test flow == 2
 
       part1, part2 = mincut_partitions(g, "A", "D", weights)
-      @test issetequal(part1, ["A", "B"])
-      @test issetequal(part2, ["C", "D"])
+      @test issetequal(part1, ["A", "B"]) || issetequal(part1, ["C", "D"])
+      @test issetequal(vcat(part1, part2), vertices(g))
 
       part1, part2 = mincut_partitions(g, weights)
-      @test issetequal(part1, ["C", "D"])
-      @test issetequal(part2, ["A", "B"])
+      @test issetequal(part1, ["A", "B"]) || issetequal(part1, ["C", "D"])
+      @test issetequal(vcat(part1, part2), vertices(g))
     end
   end
   @testset "dijkstra" begin
@@ -557,6 +556,8 @@ end
     add_edge!(g, "D" => "G")
     add_edge!(g, "E" => "G")
     t = topological_sort_by_dfs(g)
-    @test t == ["C", "B", "E", "A", "D", "G", "F"]
+    for e in edges(g)
+      @test findfirst(x -> x == src(e), t) < findfirst(x -> x == dst(e), t)
+    end
   end
 end
diff --git a/test/test_namedgraphgenerators.jl b/test/test_namedgraphgenerators.jl
index 7eb7813..6b086e7 100644
--- a/test/test_namedgraphgenerators.jl
+++ b/test/test_namedgraphgenerators.jl
@@ -30,7 +30,6 @@ using Random
   @test dims[1] > dims[2]
 
   g = triangular_lattice_graph(2, 1)
-  @show g
   dims = maximum(vertices(g))
   @test dims[1] > dims[2]
 
diff --git a/test/test_trees_and_forests.jl b/test/test_trees_and_forests.jl
new file mode 100644
index 0000000..ac84342
--- /dev/null
+++ b/test/test_trees_and_forests.jl
@@ -0,0 +1,38 @@
+using Test
+using Graphs
+using NamedGraphs
+using NamedGraphs: forest_cover, spanning_tree
+
+module TestTreesAndForests
+using NamedGraphs
+using NamedGraphs: hexagonal_lattice_graph, triangular_lattice_graph
+gs = [
+  ("Chain", named_grid((6, 1))),
+  ("Cubic Lattice", named_grid((3, 3, 3))),
+  ("Hexagonal Grid", hexagonal_lattice_graph(6, 6)),
+  ("Comb Tree", named_comb_tree((4, 4))),
+  ("Square lattice", named_grid((10, 10))),
+  ("Triangular Grid", triangular_lattice_graph(5, 5; periodic=true)),
+]
+algs = [NamedGraphs.BFS(), NamedGraphs.DFS(), NamedGraphs.RandomBFS()]
+end
+
+@testset "Test Spanning Trees $g_string, $alg" for (g_string, g) in TestTreesAndForests.gs,
+  alg in TestTreesAndForests.algs
+
+  s_tree = spanning_tree(g; alg)
+  @test is_tree(s_tree)
+  @test Set(vertices(s_tree)) == Set(vertices(g))
+  @test issubset(Set(edges(s_tree)), Set(edges(g)))
+end
+
+@testset "Test Forest Cover $g_string" for (g_string, g) in TestTreesAndForests.gs
+  cover = forest_cover(g)
+  cover_edges = reduce(vcat, edges.(cover))
+  @test issetequal(cover_edges, edges(g))
+  @test all(issetequal(vertices(forest), vertices(g)) for forest in cover)
+  for forest in cover
+    trees = NamedGraph[forest[vs] for vs in connected_components(forest)]
+    @test all(is_tree.(trees))
+  end
+end