diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
new file mode 100644
index 0000000..d1b360e
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
@@ -0,0 +1,86 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations to test well-balancedness
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 0.7,
+                                            rhos = (0.8, 0.9, 1.0))
+
+"""
+    initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations1D)
+
+Initial condition to test well balanced with a bottom topography from Fjordholm
+- Ulrik Skre Fjordholm (2012)
+  Energy conservative and stable schemes for the two-layer shallow water equations.
+  [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039)
+"""
+function initial_condition_fjordholm_well_balanced(x, t,
+                                                   equations::ShallowWaterMultiLayerEquations1D)
+    inicenter = 0.5
+    x_norm = x[1] - inicenter
+    r = abs(x_norm)
+
+    H = [0.7, 0.6, 0.5]
+    v = [0.0, 0.0, 0.0]
+    b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0
+    
+    return prim2cons(SVector(H..., v..., b), equations)
+end
+
+initial_condition = initial_condition_fjordholm_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl
index 3d74b5a..8546d9a 100644
--- a/src/equations/shallow_water_multilayer_1d.jl
+++ b/src/equations/shallow_water_multilayer_1d.jl
@@ -412,4 +412,14 @@ end
     setindex!(f, f_h, i)
     setindex!(f, f_hv, i + nlayers(equations))
 end
+
+# Calculate the error for the "lake-at-rest" test case where H = ∑h+b should
+# be a constant value over time
+@inline function Trixi.lake_at_rest_error(u,
+                                          equations::ShallowWaterMultiLayerEquations1D)
+    h = waterheight(u, equations)
+    b = u[end]
+
+    return abs(equations.H0 - (sum(h) + b))
+end
 end # @muladd
diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl
index 72f2548..16acee2 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -503,6 +503,37 @@ end # 2LSWE
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
         end
     end
+    @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_well_balanced.jl"),
+                            l2=[
+                                1.2015338985485869e-17,
+                                1.4450281975802518e-17,
+                                0.001002881985401688,
+                                5.6968770683622844e-18,
+                                5.8630384557660415e-18,
+                                1.938143907661441e-17,
+                                0.0010028819854016856,
+                            ],
+                            linf=[
+                                8.326672684688674e-17,
+                                1.3877787807814457e-16,
+                                0.005119880996670767,
+                                3.2043414465411084e-17,
+                                3.288278489632501e-17,
+                                1.6336799907854548e-16,
+                                0.005119880996670708,
+                            ],
+                            tspan=(0.0, 0.25))
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
 end # MLSWE
 end # TreeMesh1D