diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl
index 8a9384f..05af3f3 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl
@@ -6,7 +6,7 @@ using TrixiShallowWater
 ###############################################################################
 # Semidiscretization of the multilayer shallow water equations with three layers
 
-equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 10.0,
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.1,
                                               rhos = (0.9, 1.0, 1.1))
 
 initial_condition = initial_condition_convergence_test
@@ -25,7 +25,7 @@ solver = DGSEM(polydeg = 3,
 coordinates_min = 0.0
 coordinates_max = sqrt(2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 4,
+                initial_refinement_level = 2,
                 n_cells_max = 10_000,
                 periodicity = true)
 
@@ -50,12 +50,16 @@ save_solution = SaveSolutionCallback(interval = 500,
                                      save_initial_solution = true,
                                      save_final_solution = true)
 
-callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
 
-# use a Runge-Kutta method with automatic (error based) time step size control
-sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+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/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl
new file mode 100644
index 0000000..7b070fb
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl
@@ -0,0 +1,64 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with three layers
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1,
+                                              rhos = (0.9, 1.0, 1.1))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# 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, 0.0)
+coordinates_max = (sqrt(2.0), sqrt(2.0))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                n_cells_max = 100_000,
+                periodicity = true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+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/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl
new file mode 100644
index 0000000..ad61a69
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl
@@ -0,0 +1,134 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break test with a 
+# discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0,
+                                              rhos = (0.9, 0.95, 1.0))
+
+# This academic testcase sets up a discontinuous bottom topography 
+# function and initial condition to test entropy conservation. 
+
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    # Bottom topography
+    b = 0.3 * exp(-0.5 * ((x[1])^2 + (x[2])^2))
+
+    if x[1] < 0.0
+        H = SVector(1.0, 0.8, 0.6)
+    else
+        H = SVector(0.9, 0.7, 0.5)
+        b += 0.1
+    end
+
+    v1 = zero(H)
+    v2 = zero(H)
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+boundary_conditions = boundary_condition_slip_wall
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a non-periodic mesh with wall boundary conditions
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = false)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+###############################################################################
+#= 
+Workaround for TreeMesh2D to set true discontinuities for debugging and testing. 
+Essentially, this is a slight augmentation of the `compute_coefficients` where the `x` node values
+passed here are slightly perturbed in order to set a true discontinuity that avoids the doubled 
+value of the LGL nodes at a particular element interface.
+=#
+
+# Point to the data we want to augment
+u = Trixi.wrap_array(ode.u0, semi)
+# Reset the initial condition
+for element in eachelement(semi.solver, semi.cache)
+    for i in eachnode(semi.solver), j in eachnode(semi.solver)
+        x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations,
+                                       semi.solver, i, j, element)
+        # Changing the node positions passed to the initial condition by the minimum
+        # amount possible with the current type of floating point numbers allows setting
+        # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump`
+        # works if the jump location `x_jump` is at the position of an interface.
+        if i == 1 && j == 1 # bottom left corner
+            x_node = SVector(nextfloat(x_node[1]), nextfloat(x_node[2]))
+        elseif i == 1 && j == nnodes(semi.solver) # top left corner
+            x_node = SVector(nextfloat(x_node[1]), prevfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) && j == 1 # bottom right corner
+            x_node = SVector(prevfloat(x_node[1]), nextfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) && j == nnodes(semi.solver) # top right corner
+            x_node = SVector(prevfloat(x_node[1]), prevfloat(x_node[2]))
+        elseif i == 1 # left boundary
+            x_node = SVector(nextfloat(x_node[1]), x_node[2])
+        elseif j == 1 # bottom boundary
+            x_node = SVector(x_node[1], nextfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) # right boundary
+            x_node = SVector(prevfloat(x_node[1]), x_node[2])
+        elseif j == nnodes(semi.solver) # top boundary
+            x_node = SVector(x_node[1], prevfloat(x_node[2]))
+        end
+
+        u_node = initial_condition_dam_break(x_node, first(tspan), equations)
+        Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
+    end
+end
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_total,
+                                                                 energy_kinetic,
+                                                                 energy_internal))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+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/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
new file mode 100644
index 0000000..3069112
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
@@ -0,0 +1,79 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a bottom topography function 
+# to test well-balancedness
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6,
+                                              rhos = (0.7, 0.8, 0.9, 1.0))
+
+# An initial condition with constant total water height, zero velocities and a bottom topography to 
+# test well-balancedness
+function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    H = SVector(0.6, 0.55, 0.5, 0.45)
+    v1 = zero(H)
+    v2 = zero(H)
+    b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ?
+         0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] -
+                                                    0.5)^2)) + 1) : 0.0)
+
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     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/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl
new file mode 100644
index 0000000..1a9fec4
--- /dev/null
+++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl
@@ -0,0 +1,65 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a periodic
+# bottom topography function (set in the initial conditions)
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1,
+                                              rhos = (0.9, 1.0, 1.1))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 6, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# This setup is for the curved, split form convergence test on a periodic domain
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
+                           joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+# Create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+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/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl
new file mode 100644
index 0000000..5e41ef1
--- /dev/null
+++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl
@@ -0,0 +1,140 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break test with a 
+# discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0,
+                                              rhos = (0.9, 0.95, 1.0))
+
+# This test case uses a special work around to setup a truly discontinuous bottom topography 
+# function and initial condition for this academic testcase of entropy conservation. First, a 
+# dummy initial_condition_dam_break is introduced to create the semidiscretization. Then the initial
+# condition is reset with the true discontinuous values from initial_condition_discontinuous_dam_break.
+
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    # Bottom topography
+    b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2))
+
+    if x[1] < sqrt(2) / 2
+        H = SVector(1.0, 0.8, 0.6)
+    else
+        H = SVector(0.9, 0.7, 0.5)
+        b += 0.1
+    end
+
+    v1 = zero(H)
+    v2 = zero(H)
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break)
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 6, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
+                           joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = false)
+
+# Boundary conditions
+boundary_condition = Dict(:Top => boundary_condition_slip_wall,
+                          :Left => boundary_condition_slip_wall,
+                          :Right => boundary_condition_slip_wall,
+                          :Bottom => boundary_condition_slip_wall)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
+                                    solver, boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.
+
+# alternative version of the initial conditinon used to setup a truly discontinuous
+# test case and initial condition.
+# In contrast to the usual signature of initial conditions, this one get passed the
+# `element_id` explicitly. In particular, this initial conditions works as intended
+# only for the specific mesh loaded above!
+function initial_condition_discontinuous_dam_break(x, t, element_id,
+                                                   equations::ShallowWaterMultiLayerEquations2D)
+    # Bottom topography
+    b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2))
+
+    # Left side of discontinuity
+    IDs = [1, 2, 5, 6, 9, 10, 13, 14]
+    if element_id in IDs
+        H = SVector(1.0, 0.8, 0.6)
+    else # Right side of discontinuity
+        H = SVector(0.9, 0.7, 0.5)
+        b += 0.1
+    end
+
+    v1 = zero(H)
+    v2 = zero(H)
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+# point to the data we want to augment
+u = Trixi.wrap_array(ode.u0, semi)
+# reset the initial condition
+for element in eachelement(semi.solver, semi.cache)
+    for j in eachnode(semi.solver), i in eachnode(semi.solver)
+        x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations,
+                                       semi.solver, i, j, element)
+        u_node = initial_condition_discontinuous_dam_break(x_node, first(tspan), element,
+                                                           equations)
+        Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
+    end
+end
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_total,
+                                                                 energy_kinetic,
+                                                                 energy_internal))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+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/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
new file mode 100644
index 0000000..986f236
--- /dev/null
+++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl
@@ -0,0 +1,77 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a bottom topography function 
+# to test well-balancedness
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6,
+                                              rhos = (0.7, 0.8, 0.9, 1.0))
+
+# An initial condition with constant total water height, zero velocities and a bottom topography to 
+# test well-balancedness
+function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    H = SVector(0.6, 0.55, 0.5, 0.45)
+    v1 = zero(H)
+    v2 = zero(H)
+    b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ?
+         0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] -
+                                                    0.5)^2)) + 1) : 0.0)
+
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 6, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# This setup is for the curved, split form convergence test on a periodic domain
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
+                           joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+# Create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     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/TrixiShallowWater.jl b/src/TrixiShallowWater.jl
index ac74b97..c3b0147 100644
--- a/src/TrixiShallowWater.jl
+++ b/src/TrixiShallowWater.jl
@@ -20,7 +20,7 @@ include("solvers/indicators.jl")
 # Export types/functions that define the public API of TrixiShallowWater.jl
 export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
-       ShallowWaterMultiLayerEquations1D
+       ShallowWaterMultiLayerEquations1D, ShallowWaterMultiLayerEquations2D
 
 export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,
        min_max_speed_chen_noelle, flux_hll_chen_noelle,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index a36f993..89713b9 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -17,6 +17,7 @@ include("shallow_water_two_layer_2d.jl")
 abstract type AbstractShallowWaterMultiLayerEquations{NDIMS, NVARS, NLAYERS} <:
               Trixi.AbstractEquations{NDIMS, NVARS} end
 include("shallow_water_multilayer_1d.jl")
+include("shallow_water_multilayer_2d.jl")
 
 """
     eachlayer(equations::AbstractShallowWaterMultiLayerEquations)
diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl
index 328ea85..138768f 100644
--- a/src/equations/shallow_water_multilayer_1d.jl
+++ b/src/equations/shallow_water_multilayer_1d.jl
@@ -104,6 +104,7 @@ end
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations1D) = True()
+
 function Trixi.varnames(::typeof(cons2cons),
                         equations::ShallowWaterMultiLayerEquations1D)
     waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations)))
@@ -137,7 +138,7 @@ function Trixi.initial_condition_convergence_test(x, t,
     end
 
     # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
-    ω = 2.0 * pi * sqrt(2.0)
+    ω = pi * sqrt(2.0)
 
     H1 = 4.0 + 0.1 * cos(ω * x[1] + t)
     H2 = 2.0 + 0.1 * sin(ω * x[1] + t)
@@ -145,7 +146,7 @@ function Trixi.initial_condition_convergence_test(x, t,
     v1 = 0.9
     v2 = 1.0
     v3 = 1.1
-    b = 1.0 + 0.1 * cos(2.0 * ω * x[1])
+    b = 1.0 + 0.1 * cos(ω * x[1])
 
     return prim2cons(SVector(H1, H2, H3, v1, v2, v3, b), equations)
 end
@@ -162,39 +163,43 @@ in non-periodic domains).
                                                      equations::ShallowWaterMultiLayerEquations1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocity is taken to be constant
-    ω = 2 * pi * sqrt(2.0)
+    ω = pi * sqrt(2.0)
+    g = equations.gravity
 
     du1 = -0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω) +
           0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)
+
     du2 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) +
           0.1 * sin(t + x[1] * ω) * ω +
           0.1 * cos(t + x[1] * ω) * ω
+
     du3 = -0.1 * sin(t + x[1] * ω) +
-          1.1 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω)
+          1.1 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω)
+
     du4 = 0.9 * (-0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) +
           0.81 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
-          10.0 * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
+          g * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
           (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
-          (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
+          0.1 * g * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
           cos(t + x[1] * ω) * ω
+
     du5 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) +
           0.1 * sin(t + x[1] * ω) * ω +
           0.1 * cos(t + x[1] * ω) * ω +
-          10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) *
+          g * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) *
+          (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) +
+          g * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) *
           (-0.1 * sin(t + x[1] * ω) * ω +
-           0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) +
-          10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) *
-          (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω)
+           0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω))
+
     du6 = -0.11 * sin(t + x[1] * ω) +
-          1.21 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0x[1] * ω) * ω) +
-          10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
-          (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω) +
-          10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
-          (0.8181818181818181 *
-           (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
-           0.9090909090909091 *
-           (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) -
-           0.2 * sin(2.0 * x[1] * ω) * ω)
+          1.21 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) +
+          g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
+          (-0.1 * sin(x[1] * ω) * ω +
+           1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) +
+           0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) +
+          g * (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[1] * ω)) *
+          (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω)
 
     return SVector(du1, du2, du3, du4, du5, du6, zero(eltype(u)))
 end
diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl
new file mode 100644
index 0000000..da26309
--- /dev/null
+++ b/src/equations/shallow_water_multilayer_2d.jl
@@ -0,0 +1,855 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    ShallowWaterMultiLayerEquations2D(gravity, H0, rhos)
+
+Multi-Layer Shallow Water equations (MLSWE) in two space dimension. The equations are given by
+```math
+\left\{
+	\begin{aligned}			
+		&\partial_t h_m + \partial_x h_mv_m = 0,\\
+		&\partial h_mv1_m + \partial_x h_mv1_m^2 + \partial_y h_mv1_mv2_m = -gh_m\partial_x \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k<m}\frac{\rho_k}{\rho_m}h_k \bigg)
+        &\partial h_mv2_m + \partial_x h_mv1_mv2_m + \partial_y h_mv2_m^2 = -gh_m\partial_y \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k<m}\frac{\rho_k}{\rho_m}h_k \bigg)
+	\end{aligned}
+\right.
+```
+
+where ``m = 1, 2, ..., M`` is the layer index and the unknown variables are the water height ``h`` and
+the velocities ``v1, v2`` in both spatial dimensions .  Furthermore, ``g`` denotes the gravitational 
+constant, ``b(x)`` the bottom topography and ``\rho_m`` the m-th layer density, that must be chosen such that 
+``\rho_1 < \rho_2 < ... < \rho_M``, to ensure that different layers are ordered from top to bottom, with 
+increasing density.
+
+We use a specific formulation of the system, where the pressure term is reformulated as a 
+nonconservative term, which has some benefits for the design of well-balanced schemes.
+
+The additional quantity ``H_0`` is also available to store a reference value for the total water
+height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
+
+The bottom topography function ``b(x)`` is set inside the initial condition routine
+for a particular problem setup.
+
+In addition to the unknowns, Trixi currently stores the bottom topography values at the
+approximation points despite being fixed in time. This is done for convenience of computing the
+bottom topography gradients on the fly during the approximation as well as computing auxiliary
+quantities like the total water height ``H`` or the entropy variables.
+This affects the implementation and use of these equations in various ways:
+* The flux values corresponding to the bottom topography must be zero.
+* The bottom topography values must be included when defining initial conditions, boundary
+  conditions or source terms.
+* [`AnalysisCallback`](@ref) analyzes this variable.
+* Trixi's visualization tools will visualize the bottom topography by default.
+
+A good introduction for the MLSWE is available in Chapter 12 of the book:
+    - Benoit Cushman-Roisin (2011)\
+      Introduction to geophyiscal fluid dynamics: physical and numerical aspects\
+      <https://www.sciencedirect.com/bookseries/international-geophysics/vol/101/suppl/C>\
+      ISBN: 978-0-12-088759-0
+"""
+struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <:
+       AbstractShallowWaterMultiLayerEquations{2, NVARS, NLAYERS}
+    gravity::RealT   # gravitational constant
+    H0::RealT        # constant "lake-at-rest" total water height
+    rhos::SVector{NLAYERS, RealT} # Vector of layer densities
+
+    function ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity::RealT,
+                                                                      H0::RealT,
+                                                                      rhos::SVector{NLAYERS,
+                                                                                    RealT}) where {
+                                                                                                   NVARS,
+                                                                                                   NLAYERS,
+                                                                                                   RealT <:
+                                                                                                   Real
+                                                                                                   }
+        # Ensure that layer densities are all positive and in increasing order
+        issorted(rhos) ||
+            throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])"))
+        min(rhos...) > 0 || throw(ArgumentError("densities must be positive"))
+
+        new(gravity, H0, rhos)
+    end
+end
+
+# Allow for flexibility to set the gravitational constant within an elixir depending on the
+# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values.
+# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest"
+# well-balancedness test cases. 
+function ShallowWaterMultiLayerEquations2D(; gravity_constant,
+                                           H0 = zero(gravity_constant), rhos)
+
+    # Promote all variables to a common type
+    _rhos = promote(rhos...)
+    RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0))
+    __rhos = SVector(map(RealT, _rhos))
+
+    # Extract number of layers and variables
+    NLAYERS = length(rhos)
+    NVARS = 3 * NLAYERS + 1
+
+    return ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity_constant,
+                                                                    H0, __rhos)
+end
+
+@inline function Base.real(::ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}) where {
+                                                                                              NVARS,
+                                                                                              NLAYERS,
+                                                                                              RealT <:
+                                                                                              Real
+                                                                                              }
+    RealT
+end
+
+Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations2D) = True()
+
+function Trixi.varnames(::typeof(cons2cons),
+                        equations::ShallowWaterMultiLayerEquations2D)
+    waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations)))
+    momentum_1 = ntuple(n -> "h" * string(n) * "_v1", Val(nlayers(equations)))
+    momentum_2 = ntuple(n -> "h" * string(n) * "_v2", Val(nlayers(equations)))
+    return (waterheight..., momentum_1..., momentum_2..., "b")
+end
+
+# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial
+# conditions
+function Trixi.varnames(::typeof(cons2prim),
+                        equations::ShallowWaterMultiLayerEquations2D)
+    total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations)))
+    velocity_1 = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations)))
+    velocity_2 = ntuple(n -> "v" * string(n) * "_2", Val(nlayers(equations)))
+    return (total_layer_height..., velocity_1..., velocity_2..., "b")
+end
+
+# Set initial conditions at physical location `x` for time `t`
+"""
+    initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D)
+
+A smooth initial condition for a three-layer configuration used for convergence tests in combination with
+[`source_terms_convergence_test`](@ref) (and
+[`Trixi.BoundaryConditionDirichlet`](@extref) in non-periodic domains).
+"""
+function Trixi.initial_condition_convergence_test(x, t,
+                                                  equations::ShallowWaterMultiLayerEquations2D)
+    # Check if the equation has been set up with three layers
+    if nlayers(equations) != 3
+        throw(ArgumentError("This initial condition is only valid for three layers"))
+    end
+
+    # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
+    ω = pi * sqrt(2.0)
+
+    H3 = 1.5 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t)
+    H2 = 2.0 + 0.1 * sin(ω * x[1] + t) + 0.1 * sin(ω * x[2] + t)
+    H1 = 4.0 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t)
+    H = (H1, H2, H3)
+    v1 = (0.8, 0.8, 0.8)
+    v2 = (1.0, 1.0, 1.0)
+
+    b = 1.0 + 0.1 * cos(ω * x[1]) + 0.1 * cos(ω * x[2])
+
+    return prim2cons(SVector(H..., v1..., v2..., b), equations)
+end
+
+"""
+    source_terms_convergence_test(u, x, t, equations::ShallowWaterMultiLayerEquations2D)
+
+Source terms used for convergence tests with a three-layer configuration in combination with
+[`initial_condition_convergence_test`](@ref)
+(and [`Trixi.BoundaryConditionDirichlet`](@extref)
+in non-periodic domains).
+"""
+@inline function Trixi.source_terms_convergence_test(u, x, t,
+                                                     equations::ShallowWaterMultiLayerEquations2D)
+    # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
+    # this manufactured solution velocity is taken to be constant
+    ω = pi * sqrt(2.0)
+    g = equations.gravity
+
+    du1 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+           0.1 * sin(t + x[2] * ω) -
+           0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω +
+           0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) -
+           0.1 * sin(t + x[2] * ω) * ω)
+
+    du2 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) +
+           0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω +
+           0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) +
+           0.1 * sin(t + x[2] * ω) * ω)
+
+    du3 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) +
+           0.1 * sin(x[2] * ω) * ω +
+           0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) -
+           0.1 * sin(t + x[2] * ω) * ω)
+
+    du4 = (0.8 * (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+            0.1 * sin(t + x[2] * ω) -
+            0.1 * cos(t + x[1] * ω)) +
+           0.8 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) +
+           0.8^2 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
+           g *
+           (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+            0.1 * sin(t + x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
+           0.1 * g *
+           (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+            0.1 * sin(t + x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω)
+
+    du5 = (0.8 * (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) +
+            0.1 * sin(t + x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) +
+           0.8 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) +
+           0.8^2 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) +
+           g *
+           (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) +
+            0.1 * sin(t + x[2] * ω) -
+            0.1 * cos(t + x[1] * ω)) *
+           (-0.1 * sin(t + x[1] * ω) * ω +
+            0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) +
+           g *
+           (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) +
+            0.1 * sin(t + x[2] * ω) -
+            0.1 * cos(t + x[1] * ω)) *
+           (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω))
+
+    du6 = (0.8 * (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω)) +
+           0.8 * (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) +
+           0.8^2 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) +
+           g *
+           (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) +
+           g *
+           (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (-0.1 * sin(x[1] * ω) * ω +
+            0.9 / 1.1 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) +
+            1.0 / 1.1 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω)))
+
+    du7 = (-0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+           0.1 * sin(t + x[2] * ω) -
+           0.1 * cos(t + x[1] * ω) - 0.1 * cos(t + x[2] * ω) * ω +
+           0.8 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) -
+           0.1 * sin(t + x[2] * ω) * ω +
+           0.1 * g *
+           (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+            0.1 * sin(t + x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω +
+           g *
+           (2.0 + 0.1 * cos(t + x[2] * ω) - 0.1 * sin(t + x[1] * ω) -
+            0.1 * sin(t + x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω))
+
+    du8 = (0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) + 0.1 * sin(t + x[2] * ω) +
+           0.1 * cos(t + x[1] * ω) + 0.1 * cos(t + x[2] * ω) * ω +
+           0.8 * (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) +
+           0.1 * sin(t + x[2] * ω) * ω +
+           g *
+           (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) +
+            0.1 * sin(t + x[2] * ω) -
+            0.1 * cos(t + x[1] * ω)) *
+           (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) +
+           g *
+           (0.5 - 0.1 * cos(t + x[2] * ω) + 0.1 * sin(t + x[1] * ω) +
+            0.1 * sin(t + x[2] * ω) -
+            0.1 * cos(t + x[1] * ω)) *
+           (0.9 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) -
+            0.1 * sin(t + x[2] * ω) * ω))
+
+    du9 = (-0.1 * sin(t + x[1] * ω) - 0.1 * sin(t + x[2] * ω) +
+           0.1 * sin(x[2] * ω) * ω +
+           0.8 * (0.1 * sin(x[1] * ω) * ω - 0.1 * sin(t + x[1] * ω) * ω) -
+           0.1 * sin(t + x[2] * ω) * ω +
+           g *
+           (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (0.9 / 1.1 * (-0.1 * cos(t + x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω) +
+            1.0 / 1.1 * (0.1 * cos(t + x[2] * ω) * ω + 0.1 * sin(t + x[2] * ω) * ω) -
+            0.1 * sin(x[2] * ω) * ω) +
+           g *
+           (0.5 - 0.1 * cos(x[1] * ω) + 0.1 * cos(t + x[2] * ω) - 0.1 * cos(x[2] * ω) +
+            0.1 * cos(t + x[1] * ω)) *
+           (0.1 * sin(x[2] * ω) * ω - 0.1 * sin(t + x[2] * ω) * ω))
+
+    return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, zero(eltype(u)))
+end
+
+"""
+    boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
+                                 equations::ShallowWaterMultiLayerEquations2D)
+Create a boundary state by reflecting the normal velocity component and keep
+the tangential velocity component unchanged. The boundary water height is taken from
+the internal value.
+For details see Section 9.2.5 of the book:
+- Eleuterio F. Toro (2001)
+  Shock-Capturing Methods for Free-Surface Shallow Flows
+  1st edition
+  ISBN 0471987662
+"""
+@inline function Trixi.boundary_condition_slip_wall(u_inner,
+                                                    normal_direction::AbstractVector,
+                                                    x, t,
+                                                    surface_flux_function,
+                                                    equations::ShallowWaterMultiLayerEquations2D)
+    # normalize the outward pointing direction
+    normal = normal_direction / norm(normal_direction)
+
+    # Extract internal values
+    h = waterheight(u_inner, equations)
+    h_v1, h_v2 = momentum(u_inner, equations)
+    b = u_inner[end]
+
+    # compute the normal velocity
+    u_normal = normal[1] * h_v1 + normal[2] * h_v2
+
+    # create the "external" boundary solution state
+    u_boundary = SVector(h...,
+                         (h_v1 - 2.0 * u_normal * normal[1])...,
+                         (h_v2 - 2.0 * u_normal * normal[2])...,
+                         b)
+
+    # calculate the boundary flux
+    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+
+    return flux
+end
+
+"""
+    boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
+                                 surface_flux_function, equations::ShallowWaterMultiLayerEquations2D)
+
+Should be used together with [`TreeMesh`](@ref).
+"""
+@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation,
+                                                    direction, x, t,
+                                                    surface_flux_function,
+                                                    equations::ShallowWaterMultiLayerEquations2D)
+    # Extract internal values
+    h = waterheight(u_inner, equations)
+    h_v1, h_v2 = momentum(u_inner, equations)
+    b = u_inner[end]
+
+    ## get the appropriate normal vector from the orientation
+    if orientation == 1
+        u_boundary = SVector(h..., -h_v1..., h_v2..., b)
+    else # orientation == 2
+        u_boundary = SVector(h..., h_v1..., -h_v2..., b)
+    end
+
+    # Calculate boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+    end
+
+    return flux
+end
+
+# Calculate 2D advective portion of the flux for a single point
+# Note, the bottom topography has no flux
+@inline function Trixi.flux(u, orientation::Integer,
+                            equations::ShallowWaterMultiLayerEquations2D)
+
+    # Extract waterheight and momentum and compute velocities
+    h_v1, h_v2 = momentum(u, equations)
+    v1, v2 = velocity(u, equations)
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    # Calculate fluxes in each layer
+    # The momentum flux simplifies as the pressure is included in the nonconservative term
+    if orientation == 1
+        for i in eachlayer(equations)
+            f_h = h_v1[i]
+            f_hv1 = f_h * v1[i]
+            f_hv2 = f_h * v2[i]
+            setlayer!(f, f_h, f_hv1, f_hv2, i, equations)
+        end
+    else # orientation == 2
+        for i in eachlayer(equations)
+            f_h = h_v2[i]
+            f_hv1 = f_h * v1[i]
+            f_hv2 = f_h * v2[i]
+            setlayer!(f, f_h, f_hv1, f_hv2, i, equations)
+        end
+    end
+
+    return SVector(f)
+end
+
+# Calculate 2D flux for a single point in the normal direction
+# Note, this directional vector is not normalized and the bottom topography has no flux
+@inline function Trixi.flux(u, normal_direction::AbstractVector,
+                            equations::ShallowWaterMultiLayerEquations2D)
+    # Extract waterheight and momentum and compute velocities
+    h = waterheight(u, equations)
+    v1, v2 = velocity(u, equations)
+
+    v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
+    h_v_normal = h .* v_normal
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    # Calculate fluxes in each layer
+    # The momentum flux simplifies as the pressure is included in the nonconservative term
+    for i in eachlayer(equations)
+        f_h = h_v_normal[i]
+        f_hv1 = f_h * v1[i]
+        f_hv2 = f_h * v2[i]
+        setlayer!(f, f_h, f_hv1, f_hv2, i, equations)
+    end
+
+    return SVector(f)
+end
+
+"""
+    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+    flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                     normal_direction_ll::AbstractVector,
+                                     normal_direction_average::AbstractVector,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+
+Non-symmetric path-conservative two-point flux discretizing the nonconservative (source) term
+that contains the gradients of the bottom topography and waterheights from the coupling between layers
+and the nonconservative pressure formulation [`ShallowWaterMultiLayerEquations2D`](@ref).
+
+When the bottom topography is nonzero this scheme will be well-balanced when used with [`flux_ersing_etal`](@ref).
+
+In the two-layer setting this combination is equivalent to the fluxes in:
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
+"""
+@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                        orientation::Integer,
+                                                        equations::ShallowWaterMultiLayerEquations2D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    b_rr = u_rr[end]
+    b_ll = u_ll[end]
+
+    # Compute the jumps
+    h_jump = h_rr - h_ll
+    b_jump = b_rr - b_ll
+    g = equations.gravity
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    # Compute the nonconservative flux in each layer
+    # where f_hv[i] = g * h[i] * (b + ∑h[k] + ∑σ[k] * h[k])_x and σ[k] = ρ[k] / ρ[i] denotes the 
+    # density ratio of different layers
+    for i in eachlayer(equations)
+        f_hv = g * h_ll[i] * b_jump
+        for j in eachlayer(equations)
+            if j < i
+                f_hv += g * h_ll[i] *
+                        (equations.rhos[j] / equations.rhos[i] * h_jump[j])
+            else # (i<j<nlayers) nonconservative formulation of the pressure
+                f_hv += g * h_ll[i] * h_jump[j]
+            end
+        end
+
+        if orientation == 1
+            setindex!(f, f_hv, i + nlayers(equations))
+        else # orientation == 2
+            setindex!(f, f_hv, i + 2 * nlayers(equations))
+        end
+    end
+
+    return SVector(f)
+end
+
+@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                        normal_direction_ll::AbstractVector,
+                                                        normal_direction_average::AbstractVector,
+                                                        equations::ShallowWaterMultiLayerEquations2D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    b_rr = u_rr[end]
+    b_ll = u_ll[end]
+
+    # Compute the jumps
+    h_jump = h_rr - h_ll
+    b_jump = b_rr - b_ll
+    g = equations.gravity
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    for i in eachlayer(equations)
+        f_h = zero(real(equations))
+        f_hv = g * h_ll[i] * b_jump
+        for j in eachlayer(equations)
+            if j < i
+                f_hv += g * h_ll[i] *
+                        (equations.rhos[j] / equations.rhos[i] * h_jump[j])
+            else # (i<j<nlayers) nonconservative formulation of the pressure
+                f_hv += g * h_ll[i] * h_jump[j]
+            end
+        end
+        setlayer!(f, f_h, f_hv * normal_direction_average[1],
+                  f_hv * normal_direction_average[2], i, equations)
+    end
+
+    return (SVector(f))
+end
+
+"""
+    flux_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+    flux_ersing_etal(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::ShallowWaterMultiLayerEquations2D)
+
+Total energy conservative (mathematical entropy for MLSWE) split form,
+without the hydrostatic pressure.
+When the bottom topography is nonzero this scheme will be well-balanced when used with the 
+nonconservative [`flux_nonconservative_ersing_etal`](@ref).
+
+To obtain an entropy stable formulation the `surface_flux` can be set as
+`FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()), flux_nonconservative_ersing_etal`.
+
+In the two-layer setting this combination is equivalent to the fluxes in:
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2)
+"""
+@inline function flux_ersing_etal(u_ll, u_rr,
+                                  orientation::Integer,
+                                  equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_v1_ll, h_v2_ll = momentum(u_ll, equations)
+    h_v1_rr, h_v2_rr = momentum(u_rr, equations)
+
+    # Get the velocities on either side
+    v1_ll, v2_ll = velocity(u_ll, equations)
+    v1_rr, v2_rr = velocity(u_rr, equations)
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    # Calculate fluxes in each layer
+    for i in eachlayer(equations)
+        # Compute averages
+        v1_avg = 0.5 * (v1_ll[i] + v1_rr[i])
+        v2_avg = 0.5 * (v2_ll[i] + v2_rr[i])
+        h_v1_avg = 0.5 * (h_v1_ll[i] + h_v1_rr[i])
+        h_v2_avg = 0.5 * (h_v2_ll[i] + h_v2_rr[i])
+
+        # Compute fluxes
+        # The momentum flux simplifies as the pressure is included in the nonconservative term.
+        if orientation == 1
+            f_h = h_v1_avg
+            f_hv1 = f_h * v1_avg
+            f_hv2 = f_h * v2_avg
+        else # orientation == 2
+            f_h = h_v2_avg
+            f_hv1 = f_h * v1_avg
+            f_hv2 = f_h * v2_avg
+        end
+        setlayer!(f, f_h, f_hv1, f_hv2, i, equations)
+    end
+
+    return SVector(f)
+end
+
+@inline function flux_ersing_etal(u_ll, u_rr,
+                                  normal_direction::AbstractVector,
+                                  equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_v1_ll, h_v2_ll = momentum(u_ll, equations)
+    h_v1_rr, h_v2_rr = momentum(u_rr, equations)
+
+    # Get the velocities on either side
+    v1_ll, v2_ll = velocity(u_ll, equations)
+    v1_rr, v2_rr = velocity(u_rr, equations)
+
+    # Initialize flux vector
+    f = zero(MVector{3 * nlayers(equations) + 1, real(equations)})
+
+    # Calculate fluxes in each layer
+    for i in eachlayer(equations)
+        # Compute averages
+        v1_avg = 0.5 * (v1_ll[i] + v1_rr[i])
+        v2_avg = 0.5 * (v2_ll[i] + v2_rr[i])
+        h_v1_avg = 0.5 * (h_v1_ll[i] + h_v1_rr[i])
+        h_v2_avg = 0.5 * (h_v2_ll[i] + h_v2_rr[i])
+
+        # Compute fluxes
+        # The momentum flux simplifies as the pressure is included in the nonconservative term.
+        f_h = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2]
+        f_hv1 = f_h * v1_avg
+        f_hv2 = f_h * v2_avg
+
+        setlayer!(f, f_h, f_hv1, f_hv2, i, equations)
+    end
+
+    return SVector(f)
+end
+
+# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom
+# topography
+@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
+                                                              orientation_or_normal_direction,
+                                                              equations::ShallowWaterMultiLayerEquations2D)
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations)
+    diss = -0.5 * λ * (u_rr - u_ll)
+    return SVector(@views diss[1:(end - 1)]..., zero(eltype(u_ll)))
+end
+
+@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterMultiLayerEquations2D)
+    h = waterheight(u, equations)
+    h_v1, h_v2 = momentum(u, equations)
+
+    # Calculate averaged velocity of both layers
+    H = sum(h)
+    v1_m = sum(h_v1) / H
+    v2_m = sum(h_v2) / H
+    c = sqrt(equations.gravity * H)
+
+    return abs(v1_m) + c, abs(v2_m) + c
+end
+
+# Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the
+# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate
+# eigenvalues as there is no simple way to calculate them analytically.
+@inline function Trixi.max_abs_speed_naive(u_ll, u_rr,
+                                           orientation::Integer,
+                                           equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+
+    # Get the momentum quantities in the appropriate direction
+    if orientation == 1
+        h_v_ll, _ = momentum(u_ll, equations)
+        h_v_rr, _ = momentum(u_rr, equations)
+    else
+        _, h_v_ll = momentum(u_ll, equations)
+        _, h_v_rr = momentum(u_rr, equations)
+    end
+
+    # Get the averaged velocity
+    v_m_ll = sum(h_v_ll) / sum(h_ll)
+    v_m_rr = sum(h_v_rr) / sum(h_rr)
+
+    # Calculate the wave celerity on the left and right
+    c_ll = sqrt(equations.gravity * sum(h_ll))
+    c_rr = sqrt(equations.gravity * sum(h_rr))
+
+    return max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr)
+end
+
+@inline function Trixi.max_abs_speed_naive(u_ll, u_rr,
+                                           normal_direction::AbstractVector,
+                                           equations::ShallowWaterMultiLayerEquations2D)
+    # Unpack left and right state
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    h_v1_ll, h_v2_ll = momentum(u_ll, equations)
+    h_v1_rr, h_v2_rr = momentum(u_rr, equations)
+
+    # Get the averaged velocity
+    v1_m_ll = sum(h_v1_ll) / sum(h_ll)
+    v2_m_ll = sum(h_v2_ll) / sum(h_ll)
+    v1_m_rr = sum(h_v1_rr) / sum(h_rr)
+    v2_m_rr = sum(h_v2_rr) / sum(h_rr)
+
+    # Compute velocity in the normal direction
+    v_dot_n_ll = v1_m_ll * normal_direction[1] + v2_m_ll * normal_direction[2]
+    v_dot_n_rr = v1_m_rr * normal_direction[1] + v2_m_rr * normal_direction[2]
+
+    # Calculate the wave celerity on the left and right
+    c_ll = sqrt(equations.gravity * sum(h_ll))
+    c_rr = sqrt(equations.gravity * sum(h_rr))
+
+    # The normal velocities are already scaled by the norm
+    return (max(abs(v_dot_n_ll), abs(v_dot_n_rr)) +
+            max(c_ll, c_rr) * norm(normal_direction))
+end
+
+# Convert conservative variables to primitive
+@inline function Trixi.cons2prim(u, equations::ShallowWaterMultiLayerEquations2D)
+    # Extract waterheight
+    h = waterheight(u, equations)
+    b = u[end]
+
+    # Initialize total layer height
+    H = MVector{nlayers(equations), real(equations)}(undef)
+    for i in reverse(eachlayer(equations))
+        if i == nlayers(equations)
+            setindex!(H, h[i] + b, i)
+        else
+            setindex!(H, h[i] + H[i + 1], i)
+        end
+    end
+
+    v1, v2 = velocity(u, equations)
+    return SVector{3 * nlayers(equations) + 1, real(equations)}(H..., v1..., v2..., b)
+end
+
+# Convert primitive to conservative variables
+@inline function Trixi.prim2cons(prim, equations::ShallowWaterMultiLayerEquations2D)
+    # To extract the total layer height and velocity we reuse the waterheight and momentum functions 
+    # from the conservative variables.
+    H = waterheight(prim, equations)    # For primitive variables this extracts the total layer height
+    v1, v2 = momentum(prim, equations)  # For primitive variables this extracts the velocities
+    b = prim[end]
+
+    # Calculate waterheight
+    h = MVector{nlayers(equations), real(equations)}(undef)
+    for i in eachlayer(equations)
+        if i < nlayers(equations)
+            h[i] = H[i] - H[i + 1]
+        else
+            # The lowest layer is measured from the bottom topography
+            h[i] = H[i] - b
+        end
+    end
+
+    # Calculate momentum
+    h_v1 = SVector{nlayers(equations), real(equations)}(h[i] * v1[i]
+                                                        for i in eachlayer(equations))
+    h_v2 = SVector{nlayers(equations), real(equations)}(h[i] * v2[i]
+                                                        for i in eachlayer(equations))
+
+    return SVector{3 * nlayers(equations) + 1, real(equations)}(h..., h_v1..., h_v2...,
+                                                                b)
+end
+
+# Convert conservative variables to entropy variables
+# Note, only the first `3 * nlayers(equations)` are the entropy variables
+# The last entry still just carries the bottom topography values for convenience
+@inline function Trixi.cons2entropy(u, equations::ShallowWaterMultiLayerEquations2D)
+    # Extract conservative variables and compute velocity
+    h = waterheight(u, equations)
+    b = u[end]
+    v1, v2 = velocity(u, equations)
+    g = equations.gravity
+
+    # Initialize entropy variable vector
+    w = MVector{3 * nlayers(equations) + 1, real(equations)}(undef)
+
+    # Calculate entropy variables in each layer
+    for i in eachlayer(equations)
+        # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]) - 0.5 * ρ[i] * (v1[i]^2 + v2[i]^2), 
+        # where σ[k] = ρ[k] / ρ[i] denotes the density ratio of different layers
+        w1 = equations.rhos[i] * (g * b) - 0.5 * equations.rhos[i] * (v1[i]^2 + v2[i]^2)
+        for j in eachlayer(equations)
+            if j < i
+                w1 += equations.rhos[i] * g *
+                      (equations.rhos[j] / equations.rhos[i] * h[j])
+            else # i<j<nlayers
+                w1 += equations.rhos[i] * g * h[j]
+            end
+        end
+
+        w2 = equations.rhos[i] * v1[i]
+        w3 = equations.rhos[i] * v2[i]
+
+        setlayer!(w, w1, w2, w3, i, equations)
+    end
+    setindex!(w, b, nlayers(equations) * 3 + 1)
+    return SVector(w)
+end
+
+@inline function Trixi.waterheight(u, equations::ShallowWaterMultiLayerEquations2D)
+    return SVector{nlayers(equations), real(equations)}(u[i]
+                                                        for i in 1:nlayers(equations))
+end
+
+@inline function momentum(u, equations::ShallowWaterMultiLayerEquations2D)
+    h_v1 = SVector{nlayers(equations), real(equations)}(u[i]
+                                                        for i in (nlayers(equations) + 1):(2 * nlayers(equations)))
+    h_v2 = SVector{nlayers(equations), real(equations)}(u[i]
+                                                        for i in (2 * nlayers(equations) + 1):(3 * nlayers(equations)))
+    return h_v1, h_v2
+end
+
+# Helper function to extract the velocity vector from the conservative variables
+@inline function Trixi.velocity(u, equations::ShallowWaterMultiLayerEquations2D)
+    h = waterheight(u, equations)
+    h_v1, h_v2 = momentum(u, equations)
+
+    # Compute velocity
+    v1 = SVector{nlayers(equations), real(equations)}(h_v1[i] / h[i]
+                                                      for i in 1:nlayers(equations))
+    v2 = SVector{nlayers(equations), real(equations)}(h_v2[i] / h[i]
+                                                      for i in 1:nlayers(equations))
+    return v1, v2
+end
+
+# Entropy function for the multilayer shallow water equations is the total energy
+@inline function Trixi.entropy(u, equations::ShallowWaterMultiLayerEquations2D)
+    energy_total(u, equations)
+end
+
+# Calculate total energy for a conservative state `u`.
+# The total energy is composed of the kinetic energy, the hydrostatic pressure and the potential
+# energy from the bottom topography and the layer interfaces.
+@inline function Trixi.energy_total(u, equations::ShallowWaterMultiLayerEquations2D)
+    h = waterheight(u, equations)
+    v1, v2 = velocity(u, equations)
+    b = u[end]
+    g = equations.gravity
+
+    e = zero(real(equations))
+    for i in eachlayer(equations)
+        e += (equations.rhos[i] *
+              (0.5 * (h[i] * v1[i]^2 + v2[i]^2 + g * h[i]^2) + g * h[i] * b))
+        for j in 1:(i - 1)
+            e += g * equations.rhos[j] * h[j] * h[i]
+        end
+    end
+
+    return e
+end
+
+# Calculate kinetic energy for a conservative state `u`
+@inline function Trixi.energy_kinetic(u, equations::ShallowWaterMultiLayerEquations2D)
+    h = waterheight(u, equations)
+    v1, v2 = velocity(u, equations)
+
+    return (0.5 * sum(equations.rhos[i] * h[i] * (v1[i]^2 + v2[i]^2)
+                for i in eachlayer(equations)))
+end
+
+# Calculate potential energy for a conservative state `u`
+@inline function Trixi.energy_internal(u, equations::ShallowWaterMultiLayerEquations2D)
+    return energy_total(u, equations) - energy_kinetic(u, 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::ShallowWaterMultiLayerEquations2D)
+    h = waterheight(u, equations)
+    b = u[end]
+
+    return abs(equations.H0 - (sum(h) + b))
+end
+
+# Helper function to set the layer values in the flux computation
+@inline function setlayer!(f, f_h, f_hv1, f_hv2, i,
+                           equations::ShallowWaterMultiLayerEquations2D)
+    setindex!(f, f_h, i)
+    setindex!(f, f_hv1, i + nlayers(equations))
+    setindex!(f, f_hv2, i + 2 * nlayers(equations))
+    return nothing
+end
+end # @muladd
diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl
index b3c2b0e..7e6276c 100644
--- a/src/equations/shallow_water_two_layer_1d.jl
+++ b/src/equations/shallow_water_two_layer_1d.jl
@@ -84,6 +84,7 @@ function ShallowWaterTwoLayerEquations1D(; gravity_constant,
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True()
+
 function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D)
     ("h_upper", "h_v_upper",
      "h_lower", "h_v_lower", "b")
diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl
index feb7952..7b8fd6f 100644
--- a/src/equations/shallow_water_two_layer_2d.jl
+++ b/src/equations/shallow_water_two_layer_2d.jl
@@ -98,6 +98,7 @@ function ShallowWaterTwoLayerEquations2D(; gravity_constant,
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations2D) = True()
+
 function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations2D)
     ("h_upper", "h_v1_upper", "h_v2_upper", "h_lower", "h_v1_lower", "h_v2_lower", "b")
 end
diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl
index c395786..bb5e726 100644
--- a/src/equations/shallow_water_wet_dry_1d.jl
+++ b/src/equations/shallow_water_wet_dry_1d.jl
@@ -91,6 +91,7 @@ function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_con
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True()
+
 function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry1D)
     ("h", "h_v", "b")
 end
diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl
index 32c7e81..2d59294 100644
--- a/src/equations/shallow_water_wet_dry_2d.jl
+++ b/src/equations/shallow_water_wet_dry_2d.jl
@@ -94,6 +94,7 @@ function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_con
 end
 
 Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True()
+
 Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry2D) = ("h", "h_v1",
                                                                         "h_v2", "b")
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index fac1fe6..eff0178 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -83,7 +83,8 @@ isdir(outdir) && rm(outdir, recursive = true)
                                 4.4849667259339357e-14,
                                 1.9999999999999993,
                             ],
-                            tspan=(0.0, 0.25))
+                            tspan=(0.0, 0.25),
+                            atol=1e-12)
         # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl
index 2d19361..34a2585 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -145,7 +145,7 @@ isdir(outdir) && rm(outdir, recursive = true)
                             ],
                             tspan=(0.0, 0.25),
                             # Soften the tolerance as test results vary between different CPUs
-                            atol=1000 * eps())
+                            atol=2000 * eps())
         # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
@@ -480,22 +480,22 @@ end # 2LSWE
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_shallowwater_multilayer_convergence.jl"),
                             l2=[
-                                0.0024190350732389444,
-                                0.0011553525849735248,
-                                0.001427790649735925,
-                                0.0006301962973580928,
-                                0.00047816381506058254,
-                                0.0005518957180338281,
-                                0.00013912637728693124,
+                                0.0012190787287581911,
+                                0.0007900811065315603,
+                                0.00024265888138436312,
+                                0.00044695521604641624,
+                                0.0009421189161161149,
+                                0.0001681269882698457,
+                                0.000139126377287091,
                             ],
                             linf=[
-                                0.011265033642474886,
-                                0.005221927831678352,
-                                0.007355554947688914,
-                                0.0021140505291696865,
-                                0.001562655691777548,
-                                0.002409077213895827,
-                                0.0002187445586192549,
+                                0.0035625489595796367,
+                                0.0025099269565310167,
+                                0.000677952137713711,
+                                0.0019817824651813254,
+                                0.002749046992753801,
+                                0.0006393958093637853,
+                                0.00021874455861881081,
                             ],
                             tspan=(0.0, 0.25))
         # Ensure that we do not have excessive memory allocations
@@ -512,22 +512,22 @@ end # 2LSWE
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_shallowwater_multilayer_convergence.jl"),
                             l2=[
-                                0.0019005751369889307,
-                                0.0006416934635148739,
-                                0.001323748067661235,
-                                0.0002514848406032196,
-                                0.00019551914601560465,
-                                0.00029396199135642574,
-                                0.00013912637728693124,
+                                0.0007413378395106811,
+                                0.0005405760373610748,
+                                8.772497321219398e-5,
+                                0.0006969915241438401,
+                                0.0004924779641660143,
+                                0.00011816455915431224,
+                                0.000139126377287091,
                             ],
                             linf=[
-                                0.005947927726240643,
-                                0.0028032845734302647,
-                                0.003944642659854947,
-                                0.0009866355929912807,
-                                0.0009301080749861135,
-                                0.0012826136652659414,
-                                0.0002187445586192549,
+                                0.002657085141947846,
+                                0.0028102251075707296,
+                                0.00022855706065982861,
+                                0.003310637567035979,
+                                0.002490976099376596,
+                                0.0005211086991627756,
+                                0.00021874455861881081,
                             ],
                             surface_flux=(flux_lax_friedrichs,
                                           flux_nonconservative_ersing_etal),
diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl
index eec03b4..18a780a 100644
--- a/test/test_tree_2d.jl
+++ b/test/test_tree_2d.jl
@@ -442,6 +442,255 @@ end # SWE
         end
     end
 end # 2LSWE
+
+@testset "Multilayer Shallow Water" begin
+    @trixi_testset "elixir_shallowwater_multilayer_convergence.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_convergence.jl"),
+                            l2=[
+                                0.001578480912165503,
+                                0.0010725791806397446,
+                                0.00036117546296274316,
+                                0.0010107766709183063,
+                                0.0009612673392275463,
+                                0.00026390576962071974,
+                                0.001140711809558893,
+                                0.001197469124790338,
+                                0.00031384421384170294,
+                                0.00019675440964325044,
+                            ],
+                            linf=[
+                                0.006611273262461914,
+                                0.0047987596637860674,
+                                0.001476163166610922,
+                                0.004378481797734368,
+                                0.004059620398096986,
+                                0.00107253722812789,
+                                0.004650692998510841,
+                                0.005014864957408438,
+                                0.0013175223929244861,
+                                0.0004374891172380657,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with flux_lax_friedrichs" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_convergence.jl"),
+                            l2=[
+                                0.0009657284251663235,
+                                0.0007574316192008629,
+                                0.00014884847789419973,
+                                0.0008365221427534766,
+                                0.0005829573732731771,
+                                0.00013759996082561183,
+                                0.0009795167874289885,
+                                0.0007208939038902295,
+                                0.00015863343571957187,
+                                0.00019675440964325044,
+                            ],
+                            linf=[
+                                0.005343578840408814,
+                                0.005546044302341735,
+                                0.0005110705217063471,
+                                0.004897047461771331,
+                                0.004047517471484768,
+                                0.0005991894625433924,
+                                0.006191481703377466,
+                                0.005110806058533313,
+                                0.000775668423696807,
+                                0.0004374891172380657,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            tspan=(0.0, 0.25),
+                            atol=1e-11)
+        # 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
+
+    @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_well_balanced.jl"),
+                            l2=[
+                                8.476558640011521e-18,
+                                8.37932362312502e-18,
+                                1.6048527438048627e-17,
+                                0.003076923318801488,
+                                4.686259874366214e-18,
+                                4.2109008869298085e-18,
+                                3.738959845788582e-18,
+                                2.202609051659983e-17,
+                                4.6862598743662146e-18,
+                                4.2109008869298085e-18,
+                                3.738959845788581e-18,
+                                2.2026090516599833e-17,
+                                0.0030769233188014935,
+                            ],
+                            linf=[
+                                1.8041124150158794e-16,
+                                7.632783294297951e-17,
+                                2.3592239273284576e-16,
+                                0.026474051138910215,
+                                6.685363189156086e-17,
+                                5.2914444984119546e-17,
+                                6.403718148497213e-17,
+                                2.653446487129347e-16,
+                                6.685363189156086e-17,
+                                5.2914444984119546e-17,
+                                6.403718148497213e-17,
+                                2.653446487129347e-16,
+                                0.026474051138910267,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl with flux_lax_friedrichs" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_well_balanced.jl"),
+                            l2=[
+                                5.2242434650072826e-18,
+                                7.520718359606061e-18,
+                                6.0112387963921125e-18,
+                                0.0030769233188014883,
+                                5.827683385814391e-18,
+                                6.3826870710212345e-18,
+                                5.614916422687807e-18,
+                                3.6331769798446284e-17,
+                                5.8276833858143906e-18,
+                                6.382687071021235e-18,
+                                5.614916422687807e-18,
+                                3.633176979844629e-17,
+                                0.0030769233188014935,
+                            ],
+                            linf=[
+                                3.469446951953614e-17,
+                                4.85722573273506e-17,
+                                5.551115123125783e-17,
+                                0.026474051138910215,
+                                3.440924806544217e-17,
+                                4.517518549935627e-17,
+                                4.188211194330907e-17,
+                                2.034006113871492e-16,
+                                3.440924806544217e-17,
+                                4.517518549935627e-17,
+                                4.188211194330907e-17,
+                                2.034006113871492e-16,
+                                0.026474051138910267,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_dam_break.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_dam_break.jl"),
+                            l2=[
+                                0.006867196239306527,
+                                0.007212095601190101,
+                                0.014309127811972761,
+                                0.005801387118671292,
+                                0.0057555558162007536,
+                                0.006951246408270715,
+                                6.084817535430214e-5,
+                                6.068733435363927e-5,
+                                7.807791370377477e-5,
+                                0.003857583749542185,
+                            ],
+                            linf=[
+                                0.028792046975677804,
+                                0.036503344705121676,
+                                0.19267105917517297,
+                                0.028685238833612656,
+                                0.02761561086221269,
+                                0.029892140357308587,
+                                0.0004158435510054702,
+                                0.0004119057098728111,
+                                0.0007709896609887244,
+                                0.10000011323773067,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_dam_break with flux_lax_friedrichs.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_dam_break.jl"),
+                            l2=[
+                                0.007162130568658371,
+                                0.007507487998680027,
+                                0.01500636533827742,
+                                0.005593279150854283,
+                                0.005603720137202845,
+                                0.0068519457928639385,
+                                5.748086634734491e-5,
+                                5.7787474070077486e-5,
+                                7.492163249697633e-5,
+                                0.003857583749542185,
+                            ],
+                            linf=[
+                                0.037711432402825124,
+                                0.0427565001934887,
+                                0.17275386567546575,
+                                0.025348821041845597,
+                                0.02434895938442919,
+                                0.03396954814003092,
+                                0.0003026537474088657,
+                                0.00029947595638190504,
+                                0.00037146307629417057,
+                                0.10000011323773067,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            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 # TreeMesh2D
 
 # Clean up afterwards: delete TrixiShallowWater.jl output directory
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 0482e16..fa729f1 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -101,13 +101,28 @@ end
 
 @timed_testset "Multilayer shallow water conversion between conservative/entropy variables" begin
     H = (3.5, 2.5, 1.5)
-    v = (0.25, 0.1, 0.37)
+    v1 = (0.25, 0.1, 0.37)
+    v2 = (0.13, 0.2, 0.3)
     b = 0.4
 
     let equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.8,
                                                       rhos = (0.7, 0.8, 0.9))
         # test conservion between primitive and conservative variables
-        prim_vars = SVector(H..., v..., b)
+        prim_vars = SVector(H..., v1..., b)
+        cons_vars = prim2cons(prim_vars, equations)
+        @test prim_vars ≈ cons2prim(cons_vars, equations)
+
+        total_energy = energy_total(cons_vars, equations)
+        @test total_energy ≈ entropy(cons_vars, equations)
+        @test total_energy ≈
+              energy_internal(cons_vars, equations) +
+              energy_kinetic(cons_vars, equations)
+    end
+
+    let equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.8,
+                                                      rhos = (0.7, 0.8, 0.9))
+        # test conservion between primitive and conservative variables
+        prim_vars = SVector(H..., v1..., v2..., b)
         cons_vars = prim2cons(prim_vars, equations)
         @test prim_vars ≈ cons2prim(cons_vars, equations)
 
@@ -226,18 +241,33 @@ end
                                                                      0.1,
                                                                      0.2,
                                                                  ])
+    @test_throws ArgumentError ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81,
+                                                                 rhos = [
+                                                                     -1.0,
+                                                                     0.1,
+                                                                     0.2,
+                                                                 ])
     @test_throws ArgumentError ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
                                                                  rhos = [0.1, 0.3, 0.2])
+    @test_throws ArgumentError ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81,
+                                                                 rhos = [0.1, 0.3, 0.2])
     # Ensure that both tuple and array input are equivalent    
     @test ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
                                             rhos = [0.1, 0.2, 0.3]) ==
           ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
                                             rhos = (0.1, 0.2, 0.3))
+    @test ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81,
+                                            rhos = [0.1, 0.2, 0.3]) ==
+          ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81,
+                                            rhos = (0.1, 0.2, 0.3))
 
     # Check for argument error if initial condition is called with wrong number of layers
     equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
                                                   rhos = (0.1, 0.2, 0.3, 0.4))
     @test_throws ArgumentError initial_condition_convergence_test(0.0, 0.0, equations)
+    equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81,
+                                                  rhos = (0.1, 0.2, 0.3, 0.4))
+    @test_throws ArgumentError initial_condition_convergence_test(0.0, 0.0, equations)
 end
 end # Unit tests
 
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 5fcfa32..8bdd0bd 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -426,6 +426,254 @@ end # SWE
         end
     end
 end # 2LSWE
+
+@testset "Multilayer Shallow Water" begin
+    @trixi_testset "elixir_shallowwater_multilayer_convergence.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_convergence.jl"),
+                            l2=[
+                                2.845107067734722e-5,
+                                2.541332828758267e-5,
+                                1.3776844294886904e-5,
+                                3.21641609844227e-5,
+                                1.9645259635774905e-5,
+                                1.2494077534998188e-5,
+                                3.842168759551814e-5,
+                                2.3940963660476904e-5,
+                                1.5539403052020932e-5,
+                                7.244891628173627e-7,
+                            ],
+                            linf=[
+                                0.0001686597682639679,
+                                0.00016943475945224717,
+                                0.00010161333860780886,
+                                0.00018234203945666216,
+                                0.0001440335843042595,
+                                7.821014430642315e-5,
+                                0.0002159224957090089,
+                                0.00016273539529820802,
+                                9.798918563769243e-5,
+                                4.040896422807805e-6,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with flux_lax_friedrichs" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_convergence.jl"),
+                            l2=[
+                                1.918097575328871e-5,
+                                1.502019981850744e-5,
+                                5.042301180852648e-6,
+                                1.9602948809781437e-5,
+                                1.1406394406547229e-5,
+                                5.066224008439183e-6,
+                                2.3354649343891286e-5,
+                                1.433469494607502e-5,
+                                6.046742009426357e-6,
+                                7.244891628173627e-7,
+                            ],
+                            linf=[
+                                0.0001267907408812885,
+                                0.00013376989668978378,
+                                3.3611013632417475e-5,
+                                0.00013760680376173617,
+                                0.00011659018473553218,
+                                3.15196814599239e-5,
+                                0.00017067740318088553,
+                                0.00012830136581865048,
+                                3.6074428482746335e-5,
+                                4.040896422807805e-6,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_well_balanced.jl"),
+                            l2=[
+                                9.267601508692614e-18,
+                                1.3079746419456341e-17,
+                                1.1469870592306251e-17,
+                                0.0021928129262660054,
+                                5.551772799657207e-18,
+                                5.314204778937538e-18,
+                                5.555352319653062e-18,
+                                4.268241928252567e-17,
+                                6.396791646368386e-18,
+                                6.334121301133576e-18,
+                                6.427996104426942e-18,
+                                5.40279793335249e-17,
+                                0.0021928129262660024,
+                            ],
+                            linf=[
+                                2.220446049250313e-16,
+                                2.0122792321330962e-16,
+                                2.2898349882893854e-16,
+                                0.024280130945632805,
+                                6.770402623811997e-17,
+                                6.668606509370906e-17,
+                                6.930883046118016e-17,
+                                8.016754559973626e-16,
+                                1.2223596839836688e-16,
+                                1.0879445133597863e-16,
+                                8.534611214117421e-17,
+                                9.15698049490698e-16,
+                                0.024280130945632836,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl with flux_lax_friedrichs" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_well_balanced.jl"),
+                            l2=[
+                                7.747573306861602e-18,
+                                1.217768276787359e-17,
+                                9.869557273479402e-18,
+                                0.002192812926266005,
+                                5.611298116340555e-18,
+                                6.0407287687973305e-18,
+                                5.968796919323908e-18,
+                                4.954396728449881e-17,
+                                5.4141534116082716e-18,
+                                5.946309178467258e-18,
+                                6.157681863345819e-18,
+                                4.887115264719831e-17,
+                                0.0021928129262660024,
+                            ],
+                            linf=[
+                                3.469446951953614e-17,
+                                4.163336342344337e-17,
+                                1.3183898417423734e-16,
+                                0.02428013094563286,
+                                5.852343025496837e-17,
+                                6.463560282263022e-17,
+                                7.349978245944825e-17,
+                                5.309959794423629e-16,
+                                4.541940389299994e-17,
+                                6.658319696825282e-17,
+                                9.223798706573428e-17,
+                                5.425491285882212e-16,
+                                0.024280130945632836,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_dam_break.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_dam_break.jl"),
+                            l2=[
+                                0.008133096076522545,
+                                0.009059379390936938,
+                                0.0156472431836693,
+                                0.006787446774526842,
+                                0.0067495775839911685,
+                                0.007737394569507473,
+                                0.0001835435993951076,
+                                0.00018347939156198144,
+                                0.00023500246849234892,
+                                0.004003203849568449,
+                            ],
+                            linf=[
+                                0.0268742067416933,
+                                0.05243441379266703,
+                                0.1590063807340898,
+                                0.028581575095277572,
+                                0.024486565205936787,
+                                0.03953111184090555,
+                                0.002625120329243771,
+                                0.0026387938388929785,
+                                0.0038007822188701103,
+                                0.10000000026183736,
+                            ],
+                            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
+
+    @trixi_testset "elixir_shallowwater_multilayer_dam_break.jl with flux_lax_friedrichs" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_dam_break.jl"),
+                            l2=[
+                                0.0089213866827554,
+                                0.009502261925352482,
+                                0.017690193299797062,
+                                0.006845800468099525,
+                                0.006849466485572439,
+                                0.00761550757687342,
+                                0.0004491107702429648,
+                                0.00043570105077582146,
+                                0.0005832036370132536,
+                                0.004003203849568449,
+                            ],
+                            linf=[
+                                0.04276512135712615,
+                                0.054061877851283885,
+                                0.18896205511083025,
+                                0.018902398329971263,
+                                0.017901194463979406,
+                                0.025135202096835247,
+                                0.004065137704454743,
+                                0.0038290081477974037,
+                                0.005930216968236107,
+                                0.10000000026183736,
+                            ],
+                            surface_flux=(flux_lax_friedrichs,
+                                          flux_nonconservative_ersing_etal),
+                            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 # UnstructuredMesh2D
 
 # Clean up afterwards: delete TrixiShallowWater.jl output directory