diff --git a/NEWS.md b/NEWS.md
index fef7d0c8..9ba4584e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,12 @@ DispersiveShallowWater.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.10 from v0.9.x
+
+#### Changed
+
+- `convergence_test` now returns the complete convergence orders. To obtain the mean convergence rates, use `DispersiveShallowWater.calc_mean_convergence` on the convergence orders ([#285]).
+
 ## Changes when updating to v0.9 from v0.8.x
 
 #### Changed
diff --git a/src/util.jl b/src/util.jl
index db7df191..a0d5c5c5 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -124,6 +124,21 @@ function analyze_convergence(io, errors, iterations, semi::Semidiscretization, N
     return analyze_convergence(io, errors, iterations, variablenames, Ns)
 end
 
+"""
+    DispersiveShallowWater.calc_mean_convergence(eocs)
+
+Calculate the mean convergence rates from the given experimental orders of convergence `eocs`.
+The `eocs` are expected to be in the format returned by [`convergence_test`](@ref), i.e., a `Dict` where
+the keys are the error types (e.g., `:l2`, `:linf`) and the values are matrices with the EOCs for each
+variable in the columns and the iterations in the rows.
+Returns a `Dict` with the same keys as `eocs` and the mean convergence rates for all variables as values.
+"""
+function calc_mean_convergence(eocs)
+    return Dict(kind => [sum(eocs[kind][:, v]) / length(eocs[kind][:, v])
+                         for v in 1:size(eocs[kind], 2)]
+                for kind in keys(eocs))
+end
+
 # This method is called with the collected error values to actually compute and print the EOC
 function analyze_convergence(io, errors, iterations,
                              variablenames::Union{Tuple, AbstractArray}, Ns)
@@ -139,8 +154,7 @@ function analyze_convergence(io, errors, iterations,
                         log.(Ns[1:(end - 1)] ./ Ns[2:end])
                 for (kind, error) in errorsmatrix)
 
-    eoc_mean_values = Dict{Symbol, Any}()
-    eoc_mean_values[:variables] = variablenames
+    eoc_mean_values = calc_mean_convergence(eocs)
 
     for (kind, error) in errorsmatrix
         println(io, kind)
@@ -177,18 +191,15 @@ function analyze_convergence(io, errors, iterations,
         println(io, "")
 
         # Print mean EOCs
-        mean_values = zeros(nvariables)
         for v in 1:nvariables
-            mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v])
             @printf(io, "%-15s", "mean")
-            @printf(io, "%-10.2f", mean_values[v])
+            @printf(io, "%-10.2f", eoc_mean_values[kind][v])
         end
-        eoc_mean_values[kind] = mean_values
         println(io, "")
         println(io, "-"^100)
     end
 
-    return eoc_mean_values, errorsmatrix
+    return eocs, errorsmatrix
 end
 
 function extract_initial_N(example, kwargs)
diff --git a/test/test_unit.jl b/test/test_unit.jl
index eeba6aae..5d7abe9a 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -540,19 +540,21 @@ end
 
     accuracy_orders = [2, 4, 6]
     for accuracy_order in accuracy_orders
-        eoc_mean_values, _ = convergence_test(@__MODULE__, default_example(), 2, N = 256,
-                                              tspan = (0.0, 1.0),
-                                              accuracy_order = accuracy_order)
+        eocs, _ = convergence_test(@__MODULE__, default_example(), 2, N = 256,
+                                   tspan = (0.0, 1.0),
+                                   accuracy_order = accuracy_order)
+        eoc_mean_values = DispersiveShallowWater.calc_mean_convergence(eocs)
         @test isapprox(eoc_mean_values[:l2][1], accuracy_order, atol = 0.5)
         @test isapprox(eoc_mean_values[:linf][1], accuracy_order, atol = 0.5)
         @test isapprox(eoc_mean_values[:l2][2], accuracy_order, atol = 0.5)
         @test isapprox(eoc_mean_values[:linf][2], accuracy_order, atol = 0.5)
 
-        eoc_mean_values2, _ = convergence_test(@__MODULE__, default_example(), [256, 512],
-                                               tspan = (0.0, 1.0),
-                                               accuracy_order = accuracy_order)
+        eocs2, _ = convergence_test(@__MODULE__, default_example(), [256, 512],
+                                    tspan = (0.0, 1.0),
+                                    accuracy_order = accuracy_order)
+        eoc_mean_values2 = DispersiveShallowWater.calc_mean_convergence(eocs2)
         for kind in (:l2, :linf), variable in (1, 2)
-            eoc_mean_values[kind][variable] == eoc_mean_values2[kind][variable]
+            @test eoc_mean_values[kind][variable] == eoc_mean_values2[kind][variable]
         end
     end
 end