diff --git a/Project.toml b/Project.toml
index e0ef716b..69e033b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -37,7 +37,7 @@ SciMLBase = "2.11"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"
-SummationByPartsOperators = "0.5.52"
+SummationByPartsOperators = "0.5.63"
 TimerOutputs = "0.5.7"
 TrixiBase = "0.1.3"
 julia = "1.9"
diff --git a/README.md b/README.md
index 7615b6d1..0cbbeee4 100644
--- a/README.md
+++ b/README.md
@@ -10,10 +10,11 @@
 [![DOI](https://zenodo.org/badge/635090135.svg)](https://zenodo.org/doi/10.5281/zenodo.10034636)
 
 **DispersiveShallowWater.jl** is a [Julia](https://julialang.org/) package that implements structure-preserving numerical methods for dispersive shallow water models.
-To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes for two dispersive shallow water models:
+To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes for some dispersive shallow water models:
 
 * the [BBM-BBM equations with varying bottom topography](https://iopscience.iop.org/article/10.1088/1361-6544/ac3c29),
-* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924).
+* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924),
+* the [Serre-Green-Naghdi equations](https://arxiv.org/abs/2408.02665).
 
 The semidiscretizations are based on summation by parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/).
 In order to obtain fully discrete schemes, the time integration methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) are used to solve the resulting ordinary differential equations.
diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl
new file mode 100644
index 00000000..033197d7
--- /dev/null
+++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl
@@ -0,0 +1,40 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with periodic SBP operators of accuracy order 4
+accuracy_order = 4
+solver = Solver(mesh, accuracy_order)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl
new file mode 100644
index 00000000..d22196f9
--- /dev/null
+++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl
@@ -0,0 +1,41 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 128
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier pseudospectral collocation method
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))
+solver = Solver(D1, nothing)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl
new file mode 100644
index 00000000..e095c9ff
--- /dev/null
+++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 128
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier pseudospectral collocation method
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), nnodes(mesh))
+solver = Solver(D1, nothing)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+relaxation_callback = RelaxationCallback(invariant = entropy_modified)
+# Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant
+callbacks = CallbackSet(relaxation_callback, analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl
new file mode 100644
index 00000000..8bb14375
--- /dev/null
+++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl
@@ -0,0 +1,48 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
+using SparseArrays: sparse
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with periodic upwind SBP operators of accuracy order 4
+# (note that only central-type with even order are used)
+accuracy_order = 4
+D1 = upwind_operators(periodic_derivative_operator;
+                      derivative_order = 1,
+                      accuracy_order = accuracy_order - 1,
+                      xmin = xmin(mesh), xmax = xmax(mesh),
+                      N = nnodes(mesh))
+solver = Solver(D1, nothing)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)
diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl
index c963d05c..761ca562 100644
--- a/src/DispersiveShallowWater.jl
+++ b/src/DispersiveShallowWater.jl
@@ -17,14 +17,17 @@ import SciMLBase: u_modified!
 
 @reexport using StaticArrays: SVector
 using SimpleUnPack: @unpack
-using SparseArrays: sparse
-using SummationByPartsOperators: AbstractDerivativeOperator,
+using SparseArrays: sparse, issparse
+using SummationByPartsOperators: SummationByPartsOperators,
+                                 AbstractDerivativeOperator,
                                  PeriodicDerivativeOperator, PeriodicUpwindOperators,
                                  UniformPeriodicCoupledOperator,
                                  DerivativeOperator, UpwindOperators,
                                  UniformCoupledOperator,
+                                 FourierDerivativeOperator,
                                  periodic_derivative_operator,
-                                 derivative_order, integrate, mass_matrix
+                                 derivative_order, integrate, mass_matrix,
+                                 scale_by_mass_matrix!
 import SummationByPartsOperators: grid, xmin, xmax
 using TimerOutputs: TimerOutputs, print_timer, reset_timer!
 @reexport using TrixiBase: trixi_include
@@ -41,11 +44,17 @@ include("util.jl")
 
 export examples_dir, get_examples, default_example, convergence_test
 
-export BBMBBMEquations1D, BBMBBMVariableEquations1D, SvärdKalischEquations1D,
-       SvaerdKalischEquations1D
+export AbstractShallowWaterEquations,
+       BBMBBMEquations1D, BBMBBMVariableEquations1D,
+       SvärdKalischEquations1D, SvaerdKalischEquations1D,
+       SerreGreenNaghdiEquations1D
 
-export prim2prim, prim2cons, cons2prim, waterheight_total, waterheight,
-       velocity, momentum, discharge, energy_total, entropy, lake_at_rest_error,
+export prim2prim, prim2cons, cons2prim,
+       waterheight_total, waterheight,
+       velocity, momentum, discharge,
+       gravity_constant,
+       bathymetry, still_water_surface,
+       energy_total, entropy, lake_at_rest_error,
        energy_total_modified, entropy_modified
 
 export Mesh1D, xmin, xmax, nnodes
@@ -56,6 +65,8 @@ export Semidiscretization, semidiscretize, grid
 
 export boundary_condition_periodic, boundary_condition_reflecting
 
+export bathymetry_flat
+
 export initial_condition_convergence_test,
        initial_condition_manufactured, source_terms_manufactured,
        initial_condition_manufactured_reflecting, source_terms_manufactured_reflecting,
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index ddf55a6b..414de001 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -116,7 +116,7 @@ end
 
 """
     errors(analysis_callback)
-    
+
 Return the computed errors for each timestep as a named tuple.
 The shape of each entry is (nvariables, ntimesteps).
 """
@@ -330,7 +330,10 @@ function analyze(quantity, q, t, semi::Semidiscretization)
     integrate_quantity(q -> quantity(q, semi.equations), q, semi)
 end
 
-# modified entropy from Svärd-Kalisch equations need to take the whole vector `u` for every point in space
+# The modified entropy/energy of the Svärd-Kalisch and
+# Serre-Green-Naghdi equations
+# takes the whole `q` for every point in space since it requires
+# the derivative of the velocity `v_x`.
 function analyze(quantity::Union{typeof(energy_total_modified), typeof(entropy_modified)},
                  q, t, semi::Semidiscretization)
     integrate_quantity(quantity, q, semi)
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index bd604107..594607db 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMEquations1D(gravity, D, eta0 = 0.0)
+    BBMBBMEquations1D(; gravity_constant, D = 1.0, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equations are given by
 ```math
@@ -292,11 +292,4 @@ end
     return waterheight_total(q, equations) - bathymetry(q, equations)
 end
 
-@inline function energy_total(q, equations::BBMBBMEquations1D)
-    eta, v = q
-    D = still_waterdepth(q, equations)
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
-    return e
-end
-
 @inline entropy(q, equations::BBMBBMEquations1D) = energy_total(q, equations)
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index efadd002..846a2ad5 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMVariableEquations1D(gravity, eta0 = 1.0)
+    BBMBBMVariableEquations1D(; gravity_constant, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension with spatially varying bathymetry. The equations are given by
 ```math
@@ -378,12 +378,6 @@ end
     return waterheight_total(q, equations) - bathymetry(q, equations)
 end
 
-@inline function energy_total(q, equations::BBMBBMVariableEquations1D)
-    eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
-    return e
-end
-
 @inline entropy(q, equations::BBMBBMVariableEquations1D) = energy_total(q, equations)
 
 # Calculate the error for the "lake-at-rest" test case where eta should
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index ba828f92..eba21306 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -4,6 +4,8 @@
 An abstract supertype of specific equations such as the BBM-BBM equations.
 The type parameters encode the number of spatial dimensions (`NDIMS`) and the
 number of primary variables (`NVARS`) of the physics model.
+
+See also [`AbstractShallowWaterEquations`](@ref).
 """
 abstract type AbstractEquations{NDIMS, NVARS} end
 
@@ -40,6 +42,26 @@ Common choices of the `conversion_function` are [`prim2prim`](@ref) and
 """
 function varnames end
 
+"""
+    AbstractShallowWaterEquations{NDIMS, NVARS}
+
+An abstract supertype of all equation system that contain the classical
+shallow water equations as a subsystem, e.g., the
+[`BBMBBMEquations1D`](@ref), the [`SvaerdKalischEquations1D`](@ref),
+and the [`SerreGreenNaghdiEquations1D`](@ref).
+In 1D, the shallow water equations with flat bathymetry are given by
+```math
+\begin{aligned}
+  h_t + (h v)_x &= 0,\\
+  h v_t + \frac{1}{2} g (h^2)_x + \frac{1}{2} h (v^2)_x &= 0,
+\end{aligned}
+```
+where ``h`` is the [`waterheight`](@ref),
+``v`` the [`velocity`](@ref), and
+``g`` the [`gravity_constant`](@ref).
+"""
+abstract type AbstractShallowWaterEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+
 """
     prim2prim(q, equations)
 
@@ -74,7 +96,8 @@ function cons2prim end
     waterheight_total(q, equations)
 
 Return the total waterheight of the primitive variables `q` for a given set of
-`equations`, i.e. the waterheight plus the bathymetry.
+`equations`, i.e., the [`waterheight`](@ref) plus the
+[`bathymetry`](@ref).
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
@@ -87,7 +110,7 @@ varnames(::typeof(waterheight_total), equations) = ("η",)
     waterheight(q, equations)
 
 Return the waterheight of the primitive variables `q` for a given set of
-`equations`, i.e. the waterheight above the bathymetry.
+`equations`, i.e., the waterheight above the bathymetry.
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
@@ -112,8 +135,8 @@ varnames(::typeof(velocity), equations) = ("v",)
 """
     momentum(q, equations)
 
-Return the momentum of the primitive variables `q` for a given set of
-`equations`, i.e. the waterheight times the velocity.
+Return the momentum/discharge of the primitive variables `q` for a given set of
+`equations`, i.e., the [`waterheight`](@ref) times the [`velocity`](@ref).
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
@@ -133,38 +156,129 @@ See [`momentum`](@ref).
 
 varnames(::typeof(discharge), equations) = ("P",)
 
-@inline function still_waterdepth(q, equations::AbstractEquations)
+"""
+    still_water_surface(q, equations::AbstractShallowWaterEquations)
+
+Return the still water surface ``\\eta_0`` (lake at rest)
+for a given set of `equations`.
+"""
+@inline function still_water_surface(q, equations::AbstractShallowWaterEquations)
+    return equations.eta0
+end
+
+@inline function still_waterdepth(q, equations::AbstractShallowWaterEquations)
     b = bathymetry(q, equations)
-    D = equations.eta0 - b
+    eta0 = still_water_surface(q, equations)
+    D = eta0 - b
     return D
 end
 
+"""
+    bathymetry(q, equations)
+
+Return the bathymetry of the primitive variables `q` for a given set of
+`equations`.
+
+`q` is a vector of the primitive variables at a single node, i.e., a vector
+of the correct length `nvariables(equations)`.
+"""
+function bathymetry end
+
 """
     entropy(q, equations)
 
 Return the entropy of the primitive variables `q` for a given set of
-`equations`.
+`equations`. For all [`AbstractShallowWaterEquations`](@ref), the `entropy`
+is just the [`energy_total`](@ref).
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
 """
 function entropy end
 
+function entropy(q, equations::AbstractShallowWaterEquations)
+    return energy_total(q, equations)
+end
+
 varnames(::typeof(entropy), equations) = ("U",)
 
+"""
+    gravity_constant(equations::AbstractShallowWaterEquations)
+
+Return the gravity constant ``g`` for a given set of `equations`.
+See also [`AbstractShallowWaterEquations`](@ref).
+"""
+@inline function gravity_constant(equations::AbstractShallowWaterEquations)
+    return equations.gravity
+end
+
 """
     energy_total(q, equations)
 
 Return the total energy of the primitive variables `q` for a given set of
-`equations`.
+`equations`. For all [`AbstractShallowWaterEquations`](@ref), the total
+energy is given by the sum of the kinetic and potential energy of the
+shallow water subsystem, i.e.,
+```math
+\\frac{1}{2} h v^2 + \\frac{1}{2} g \\eta^2
+```
+in 1D, where ``h`` is the [`waterheight`](@ref),
+``\\eta = h + b`` the [`waterheight_total`](@ref),
+``v`` the [`velocity`](@ref), and ``g`` the [`gravity_constant`](@ref).
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
 """
-function energy_total end
+@inline function energy_total(q, equations::AbstractShallowWaterEquations)
+    h = waterheight(q, equations)
+    eta = waterheight_total(q, equations)
+    v = velocity(q, equations)
+    return 0.5f0 * h * v^2 + 0.5f0 * gravity_constant(equations) * eta^2
+end
 
 varnames(::typeof(energy_total), equations) = ("e_total",)
 
+# The modified entropy/total energy takes the whole `q_global` for every point in space
+"""
+    energy_total_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+
+Return the modified total energy of the primitive variables `q_global` for the
+`equations`. This modified total energy is a conserved quantity and can
+contain additional terms compared to the usual [`energy_total`](@ref).
+For example, for the [`SvaerdKalischEquations1D`](@ref) and the
+[`SerreGreenNaghdiEquations1D`](@ref), it contains additional terms
+depending on the derivative of the velocity ``v_x`` modeling non-hydrostatic
+contributions.
+
+`q_global` is a vector of the primitive variables at ALL nodes.
+`cache` needs to hold the SBP operators used by the `solver` if non-hydrostatic
+terms are present.
+"""
+function energy_total_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+    # `q_global` is an `ArrayPartition` of the primitive variables at all nodes
+    @assert nvariables(equations) == length(q_global.x)
+
+    e = similar(q_global.x[begin])
+    for i in eachindex(q_global.x[begin])
+        e[i] = energy_total(get_node_vars(q_global, equations, i), equations)
+    end
+
+    return e
+end
+
+varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
+
+"""
+    entropy_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+
+Alias for [`energy_total_modified`](@ref).
+"""
+@inline function entropy_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+    energy_total_modified(q_global, equations, cache)
+end
+
+varnames(::typeof(entropy_modified), equations) = ("U_modified",)
+
 # Add methods to show some information on systems of equations.
 function Base.show(io::IO, equations::AbstractEquations)
     # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.
@@ -210,12 +324,27 @@ Default analysis integrals used by the [`AnalysisCallback`](@ref).
 """
 default_analysis_integrals(::AbstractEquations) = Symbol[]
 
+abstract type AbstractBathymetry end
+struct BathymetryFlat <: AbstractBathymetry end
+"""
+    bathymetry_flat = DispersiveShallowWater.BathymetryFlat()
+
+A singleton struct indicating a flat bathymetry.
+"""
+const bathymetry_flat = BathymetryFlat()
+
 # BBM-BBM equations
-abstract type AbstractBBMBBMEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+abstract type AbstractBBMBBMEquations{NDIMS, NVARS} <:
+              AbstractShallowWaterEquations{NDIMS, NVARS} end
 include("bbm_bbm_1d.jl")
 include("bbm_bbm_variable_bathymetry_1d.jl")
 
 # Svärd-Kalisch equations
 abstract type AbstractSvaerdKalischEquations{NDIMS, NVARS} <:
-              AbstractEquations{NDIMS, NVARS} end
+              AbstractShallowWaterEquations{NDIMS, NVARS} end
 include("svaerd_kalisch_1d.jl")
+
+# Serre-Green-Naghdi equations
+abstract type AbstractSerreGreenNaghdiEquations{NDIMS, NVARS} <:
+              AbstractShallowWaterEquations{NDIMS, NVARS} end
+include("serre_green_naghdi_flat_1d.jl")
diff --git a/src/equations/serre_green_naghdi_flat_1d.jl b/src/equations/serre_green_naghdi_flat_1d.jl
new file mode 100644
index 00000000..0858679f
--- /dev/null
+++ b/src/equations/serre_green_naghdi_flat_1d.jl
@@ -0,0 +1,456 @@
+@doc raw"""
+    SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0)
+
+Serre-Green-Naghdi system in one spatial dimension.
+The equations for flat bathymetry are given by
+```math
+\begin{aligned}
+  h_t + (h v)_x &= 0,\\
+  h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} g (h^2)_x + \frac{1}{2} h (v^2)_x + p_x &= 0,\\
+  p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx}.
+\end{aligned}
+```
+The unknown quantities of the Serre-Green-Naghdi equations are the
+total water height ``\eta = h + b`` and the velocity ``v``.
+The gravitational constant is denoted by `g` and the bottom topography
+(bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry
+is therefore given by ``h = \eta - \eta_0 + D``. The Serre-Green-Naghdi
+equations are only implemented for ``\eta_0 = 0``.
+The total water height is therefore given by ``\eta = h + b``.
+
+References for the Serre-Green-Naghdi system can be found in
+- Serre (1953)
+  Contribution â l'étude des écoulements permanents et variables dans les canaux
+  [DOI: 10.1051/lhb/1953034](https://doi.org/10.1051/lhb/1953034)
+- Green and Naghdi (1976)
+  A derivation of equations for wave propagation in water of variable depth
+  [DOI: 10.1017/S0022112076002425](https://doi.org/10.1017/S0022112076002425)
+
+The semidiscretization implemented here conserves
+- the total water mass (integral of h) as a linear invariant
+- the total momentum (integral of h v) as a nonlinear invariant
+- the total energy
+
+for periodic boundary conditions, see
+- Hendrik Ranocha and Mario Ricchiuto (2024)
+  Structure-preserving approximations of the Serre-Green-Naghdi
+  equations in standard and hyperbolic form
+  [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665)
+"""
+struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
+       AbstractSerreGreenNaghdiEquations{1, 3}
+    bathymetry::Bathymetry # type of bathymetry
+    gravity::RealT # gravitational constant
+    eta0::RealT # constant still-water surface
+end
+
+function SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0)
+    eta0 == 0.0 ||
+        @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations"
+    SerreGreenNaghdiEquations1D(bathymetry_flat, gravity_constant, eta0)
+end
+
+varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D")
+varnames(::typeof(prim2cons), ::SerreGreenNaghdiEquations1D) = ("h", "hv", "b")
+
+"""
+    initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D, mesh)
+
+A soliton solution used for convergence tests in a periodic domain.
+"""
+function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEquations1D,
+                                            mesh)
+    g = equations.gravity
+
+    # setup parameters data
+    h1 = 1.0
+    h2 = 1.2
+    c = sqrt(g * h2)
+
+    x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh)
+
+    h = h1 + (h2 - h1) * sech(x_t / 2 * sqrt(3 * (h2 - h1) / (h1^2 * h2)))^2
+    v = c * (1 - h1 / h)
+
+    return SVector(h, v, 0)
+end
+
+function create_cache(mesh,
+                      equations::SerreGreenNaghdiEquations1D{BathymetryFlat},
+                      solver,
+                      initial_condition,
+                      ::BoundaryConditionPeriodic,
+                      RealT, uEltype)
+    D = solver.D1
+
+    # create temporary storage
+    h = ones(RealT, nnodes(mesh))
+    h_x = zero(h)
+    v_x = zero(h)
+    h2_x = zero(h)
+    hv_x = zero(h)
+    v2_x = zero(h)
+    h2_v_vx_x = zero(h)
+    h_vx_x = zero(h)
+    p_x = zero(h)
+    tmp = zero(h)
+    M_h = zero(h)
+    M_h3_3 = zero(h)
+
+    if D isa PeriodicUpwindOperators
+        v_x_upwind = zero(h)
+
+        Dmat_minus = sparse(D.minus)
+
+        # Floating point errors accumulate a bit and the system matrix
+        # is not necessarily perfectly symmetric but only up to
+        # round-off errors. We wrap it here to avoid issues with the
+        # factorization.
+        @. M_h = h
+        scale_by_mass_matrix!(M_h, D)
+        @. M_h3_3 = (1 / 3) * h^3
+        scale_by_mass_matrix!(M_h3_3, D)
+        system_matrix = Symmetric(Diagonal(M_h)
+                                  +
+                                  Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus)
+        factorization = cholesky(system_matrix)
+
+        cache = (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x,
+                 h2_v_vx_x, h_vx_x, p_x, tmp,
+                 M_h, M_h3_3,
+                 D, Dmat_minus, factorization)
+    else
+        if D isa FourierDerivativeOperator
+            Dmat = Matrix(D)
+
+            cache = (; h_x, v_x, h2_x, hv_x, v2_x,
+                     h2_v_vx_x, h_vx_x, p_x, tmp,
+                     M_h, M_h3_3,
+                     D, Dmat)
+        else
+            Dmat = sparse(D)
+
+            # Floating point errors accumulate a bit and the system matrix
+            # is not necessarily perfectly symmetric but only up to
+            # round-off errors. We wrap it here to avoid issues with the
+            # factorization.
+            @. M_h = h
+            scale_by_mass_matrix!(M_h, D)
+            @. M_h3_3 = (1 / 3) * h^3
+            scale_by_mass_matrix!(M_h3_3, D)
+            system_matrix = Symmetric(Diagonal(M_h)
+                                      +
+                                      Dmat' * Diagonal(M_h3_3) * Dmat)
+            factorization = cholesky(system_matrix)
+
+            cache = (; h_x, v_x, h2_x, hv_x, v2_x,
+                     h2_v_vx_x, h_vx_x, p_x, tmp,
+                     M_h, M_h3_3,
+                     D, Dmat, factorization)
+        end
+    end
+
+    return cache
+end
+
+# Discretization that conserves
+# - the total water mass (integral of h) as a linear invariant
+# - the total momentum (integral of h v) as a nonlinear invariant
+# - the total energy
+# for periodic boundary conditions, see
+# - Hendrik Ranocha and Mario Ricchiuto (2024)
+#   Structure-preserving approximations of the Serre-Green-Naghdi
+#   equations in standard and hyperbolic form
+#   [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665)
+# TODO: Implement source terms
+# TODO: Implement variable bathymetry
+function rhs!(dq, q, t, mesh,
+              equations::SerreGreenNaghdiEquations1D,
+              initial_condition,
+              ::BoundaryConditionPeriodic,
+              source_terms::Nothing,
+              solver, cache)
+    if cache.D isa PeriodicUpwindOperators
+        rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry)
+    else
+        rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry)
+    end
+
+    return nothing
+end
+
+function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat)
+    # Unpack physical parameters and SBP operator `D` as well as the
+    # SBP operator in sparse matrix form `Dmat`
+    g = gravity_constant(equations)
+    (; D, Dmat) = cache
+
+    # `q` and `dq` are `ArrayPartition`s. They collect the individual
+    # arrays for the water height `h` and the velocity `v`.
+    h, v = q.x
+    dh, dv = dq.x
+
+    @trixi_timeit timer() "hyperbolic terms" begin
+        # Compute all derivatives required below
+        (; h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x,
+        h_vx_x, p_x, tmp, M_h, M_h3_3) = cache
+
+        mul!(h_x, D, h)
+        mul!(v_x, D, v)
+        @. tmp = h^2
+        mul!(h2_x, D, tmp)
+        @. tmp = h * v
+        mul!(hv_x, D, tmp)
+        @. tmp = v^2
+        mul!(v2_x, D, tmp)
+
+        @. tmp = h^2 * v * v_x
+        mul!(h2_v_vx_x, D, tmp)
+        @. tmp = h * v_x
+        mul!(h_vx_x, D, tmp)
+        inv6 = 1 / 6
+        @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x
+                  -
+                  inv6 * h * h2_v_vx_x
+                  -
+                  inv6 * h^2 * v * h_vx_x)
+        mul!(p_x, D, tmp)
+
+        # Plain: h_t + (h v)_x = 0
+        #
+        # Split form for energy conservation:
+        # h_t + h_x v + h v_x = 0
+        @. dh = -(h_x * v + h * v_x)
+
+        # Plain: h v_t + ... = 0
+        #
+        # Split form for energy conservation:
+        @. tmp = -(g * h2_x - g * h * h_x
+                   +
+                   0.5 * h * v2_x
+                   -
+                   0.5 * v^2 * h_x
+                   +
+                   0.5 * hv_x * v
+                   -
+                   0.5 * h * v * v_x
+                   +
+                   p_x)
+    end
+
+    @trixi_timeit timer() "assembling elliptic operator" begin
+        # The code below is equivalent to
+        #   dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp
+        # but faster since the symbolic factorization is reused.
+        # Floating point errors accumulate a bit and the system matrix
+        # is not necessarily perfectly symmetric but only up to round-off
+        # errors. We wrap it here to avoid issues with the factorization.
+        @. M_h = h
+        scale_by_mass_matrix!(M_h, D)
+        inv3 = 1 / 3
+        @. M_h3_3 = inv3 * h^3
+        scale_by_mass_matrix!(M_h3_3, D)
+        system_matrix = Symmetric(Diagonal(M_h)
+                                  +
+                                  Dmat' * Diagonal(M_h3_3) * Dmat)
+    end
+
+    @trixi_timeit timer() "solving elliptic system" begin
+        if issparse(system_matrix)
+            (; factorization) = cache
+            cholesky!(factorization, system_matrix; check = false)
+            if issuccess(factorization)
+                scale_by_mass_matrix!(tmp, D)
+                dv .= factorization \ tmp
+            else
+                # The factorization may fail if the time step is too large
+                # and h becomes negative.
+                fill!(dv, NaN)
+            end
+        else
+            factorization = cholesky!(system_matrix)
+            scale_by_mass_matrix!(tmp, D)
+            ldiv!(dv, factorization, tmp)
+        end
+    end
+
+    return nothing
+end
+
+function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat)
+    # Unpack physical parameters and SBP operator `D` as well as the
+    # SBP upwind operator in sparse matrix form `Dmat_minus`
+    g = gravity_constant(equations)
+    (; Dmat_minus) = cache
+    D_upwind = cache.D
+    D = D_upwind.central
+
+    # `q` and `dq` are `ArrayPartition`s. They collect the individual
+    # arrays for the water height `h` and the velocity `v`.
+    h, v = q.x
+    dh, dv = dq.x
+
+    @trixi_timeit timer() "hyperbolic terms" begin
+        # Compute all derivatives required below
+        (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x,
+        h2_v_vx_x, h_vx_x, p_x, tmp,
+        M_h, M_h3_3) = cache
+
+        mul!(h_x, D, h)
+        mul!(v_x, D, v)
+        mul!(v_x_upwind, D_upwind.minus, v)
+        @. tmp = h^2
+        mul!(h2_x, D, tmp)
+        @. tmp = h * v
+        mul!(hv_x, D, tmp)
+        @. tmp = v^2
+        mul!(v2_x, D, tmp)
+
+        @. tmp = h^2 * v * v_x
+        mul!(h2_v_vx_x, D, tmp)
+        @. tmp = h * v_x
+        mul!(h_vx_x, D, tmp)
+        # p_+
+        @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind
+        mul!(p_x, D_upwind.plus, tmp)
+        # p_0
+        minv6 = -1 / 6
+        @. tmp = minv6 * (h * h2_v_vx_x
+                          +
+                          h^2 * v * h_vx_x)
+        mul!(p_x, D, tmp, 1.0, 1.0)
+
+        # Plain: h_t + (h v)_x = 0
+        #
+        # Split form for energy conservation:
+        # h_t + h_x v + h v_x = 0
+        @. dh = -(h_x * v + h * v_x)
+
+        # Plain: h v_t + ... = 0
+        #
+        # Split form for energy conservation:
+        @. tmp = -(g * h2_x - g * h * h_x
+                   +
+                   0.5 * h * v2_x
+                   -
+                   0.5 * v^2 * h_x
+                   +
+                   0.5 * hv_x * v
+                   -
+                   0.5 * h * v * v_x
+                   +
+                   p_x)
+    end
+
+    @trixi_timeit timer() "assembling elliptic operator" begin
+        # The code below is equivalent to
+        #   dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp
+        # but faster since the symbolic factorization is reused.
+        # Floating point errors accumulate a bit and the system matrix
+        # is not necessarily perfectly symmetric but only up to round-off errors.
+        # We wrap it here to avoid issues with the factorization.
+        @. M_h = h
+        scale_by_mass_matrix!(M_h, D)
+        inv3 = 1 / 3
+        @. M_h3_3 = inv3 * h^3
+        scale_by_mass_matrix!(M_h3_3, D)
+        system_matrix = Symmetric(Diagonal(M_h)
+                                  +
+                                  Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus)
+    end
+
+    @trixi_timeit timer() "solving elliptic system" begin
+        (; factorization) = cache
+        cholesky!(factorization, system_matrix; check = false)
+        if issuccess(factorization)
+            scale_by_mass_matrix!(tmp, D)
+            dv .= factorization \ tmp
+        else
+            # The factorization may fail if the time step is too large
+            # and h becomes negative.
+            fill!(dv, NaN)
+        end
+    end
+
+    return nothing
+end
+
+@inline function prim2cons(q, equations::SerreGreenNaghdiEquations1D)
+    h = waterheight(q, equations)
+    v = velocity(q, equations)
+    b = bathymetry(q, equations)
+
+    hv = h * v
+    return SVector(h, hv, b)
+end
+
+@inline function cons2prim(u, equations::SerreGreenNaghdiEquations1D)
+    h, hv, b = u
+
+    eta = h + b
+    v = hv / h
+    D = equations.eta0 - b
+    return SVector(eta, v, D)
+end
+
+@inline function waterheight_total(q, equations::SerreGreenNaghdiEquations1D)
+    return q[1]
+end
+
+@inline function velocity(q, equations::SerreGreenNaghdiEquations1D)
+    return q[2]
+end
+
+@inline function bathymetry(q, equations::SerreGreenNaghdiEquations1D)
+    D = q[3]
+    return equations.eta0 - D
+end
+
+@inline function waterheight(q, equations::SerreGreenNaghdiEquations1D)
+    return waterheight_total(q, equations) - bathymetry(q, equations)
+end
+
+# The entropy/energy takes the whole `q` for every point in space
+"""
+    energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache)
+
+Return the modified total energy of the primitive variables `q_global` for the
+[`SerreGreenNaghdiEquations1D`](@ref).
+It contains an additional term containing a
+derivative compared to the usual [`energy_total`](@ref) modeling
+non-hydrostatic contributions. The [`energy_total_modified`](@ref)
+is a conserved quantity (for periodic boundary conditions).
+
+For a [`bathymetry_flat`](@ref) the total energy is given by
+```math
+\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2.
+```
+
+`q_global` is a vector of the primitive variables at ALL nodes.
+`cache` needs to hold the SBP operators used by the `solver`.
+"""
+function energy_total_modified(q_global,
+                               equations::SerreGreenNaghdiEquations1D,
+                               cache)
+    # unpack physical parameters and SBP operator `D`
+    g = equations.gravity
+    (; D, v_x) = cache
+
+    # `q_global` is an `ArrayPartition`. It collects the individual arrays for
+    # the water height `h` and the velocity `v`.
+    h, v = q_global.x
+
+    N = length(v)
+    e = zeros(eltype(q_global), N)
+
+    # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2
+    if D isa PeriodicUpwindOperators
+        mul!(v_x, D.minus, v)
+    else
+        mul!(v_x, D, v)
+    end
+
+    @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2
+
+    return e
+end
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index ae20581b..f07dc6fd 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -1,5 +1,8 @@
 @doc raw"""
-    SvaerdKalischEquations1D(gravity, eta0 = 1.0, alpha = 0.0, beta = 0.2308939393939394, gamma = 0.04034343434343434)
+    SvaerdKalischEquations1D(; gravity_constant, eta0 = 0.0,
+                               alpha = 0.0,
+                               beta = 0.2308939393939394,
+                               gamma = 0.04034343434343434)
 
 Dispersive system by Svärd and Kalisch in one spatial dimension with spatially varying bathymetry. The equations are given in conservative variables by
 ```math
@@ -327,34 +330,26 @@ end
     return waterheight_total(q, equations) - bathymetry(q, equations)
 end
 
-@inline function energy_total(q, equations::SvaerdKalischEquations1D)
-    eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
-    return e
-end
-
 @inline entropy(u, equations::SvaerdKalischEquations1D) = energy_total(u, equations)
 
 # The modified entropy/total energy takes the whole `q` for every point in space
 """
-    energy_total_modified(q, equations::SvaerdKalischEquations1D, cache)
+    energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)
 
-Return the modified total energy of the primitive variables `q` for the
+Return the modified total energy of the primitive variables `q_global` for the
 `SvaerdKalischEquations1D`. It contains an additional term containing a
 derivative compared to the usual `energy_total`. The `energy_total_modified`
 is a conserved quantity of the Svärd-Kalisch equations.
 
-`q` is a vector of the primitive variables at ALL nodes, i.e., a matrix
-of the correct length `nvariables(equations)` as first dimension and the
-number of nodes as length of the second dimension.
+`q_global` is a vector of the primitive variables at ALL nodes.
 `cache` needs to hold the first-derivative SBP operator `D1`.
 """
-@inline function energy_total_modified(q, equations::SvaerdKalischEquations1D, cache)
+@inline function energy_total_modified(q_global, equations::SvaerdKalischEquations1D, cache)
     # Need to compute new beta_hat, do not use the old one from the `cache`
-    v = q.x[2]
-    D = q.x[3]
+    v = q_global.x[2]
+    D = q_global.x[3]
     N = length(v)
-    e_modified = zeros(eltype(q), N)
+    e_modified = zeros(eltype(q_global), N)
     beta_hat = equations.beta * D .^ 3
     if cache.D1 isa PeriodicDerivativeOperator ||
        cache.D1 isa UniformPeriodicCoupledOperator
@@ -365,24 +360,12 @@ number of nodes as length of the second dimension.
         @error "unknown type of first-derivative operator: $(typeof(cache.D1))"
     end
     for i in 1:N
-        e_modified[i] = energy_total(get_node_vars(q, equations, i), equations) + tmp[i]
+        e_modified[i] = energy_total(get_node_vars(q_global, equations, i), equations) +
+                        tmp[i]
     end
     return e_modified
 end
 
-varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
-
-"""
-    entropy_modified(q, equations::SvaerdKalischEquations1D, cache)
-
-Alias for [`energy_total_modified`](@ref).
-"""
-@inline function entropy_modified(q, equations::SvaerdKalischEquations1D, cache)
-    energy_total_modified(q, equations, cache)
-end
-
-varnames(::typeof(entropy_modified), equations) = ("U_modified",)
-
 # Calculate the error for the "lake-at-rest" test case where eta should
 # be a constant value over time
 @inline function lake_at_rest_error(u, equations::SvaerdKalischEquations1D)
diff --git a/src/semidiscretization.jl b/src/semidiscretization.jl
index b13599e7..6e13dde7 100644
--- a/src/semidiscretization.jl
+++ b/src/semidiscretization.jl
@@ -144,7 +144,9 @@ function integrate_quantity!(quantity, func, q, semi::Semidiscretization)
     integrate(quantity, semi)
 end
 
-# modified entropy from Svärd-Kalisch equations need to take the whole vector `q` for every point in space
+# The entropy/energy of the Svärd-Kalisch and Serre-Green-Naghdi equations
+# takes the whole `q` for every point in space since it requires
+# the derivative of the velocity `v_x`.
 function integrate_quantity(func::Union{typeof(energy_total_modified),
                                         typeof(entropy_modified)}, q,
                             semi::Semidiscretization)
diff --git a/src/solver.jl b/src/solver.jl
index 0b9bd9e4..a510b884 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -10,9 +10,10 @@ abstract type AbstractSolver end
 
 A `struct` that holds the summation by parts (SBP) operators that are used for the spatial discretization.
 """
-struct Solver{RealT <: Real, FirstDerivative <: AbstractDerivativeOperator{RealT},
+struct Solver{RealT <: Real,
+              FirstDerivative <: AbstractDerivativeOperator{RealT},
               SecondDerivative <:
-              Union{AbstractDerivativeOperator{RealT}, AbstractMatrix{RealT}}} <:
+              Union{AbstractDerivativeOperator{RealT}, AbstractMatrix{RealT}, Nothing}} <:
        AbstractSolver
     D1::FirstDerivative
     D2::SecondDerivative
@@ -24,7 +25,8 @@ struct Solver{RealT <: Real, FirstDerivative <: AbstractDerivativeOperator{RealT
                                                                                            SecondDerivative
                                                                                            }
         @assert derivative_order(D1) == 1
-        if D2 isa AbstractDerivativeOperator
+        if D2 isa AbstractDerivativeOperator &&
+           !(D2 isa SummationByPartsOperators.FourierPolynomialDerivativeOperator)
             @assert derivative_order(D2) == 2
         end
         new(D1, D2)
@@ -48,14 +50,16 @@ end
 """
     Solver(D1, D2)
 
-Create a solver, where `D1` is an `AbstractDerivativeOperator`  from [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
-of first `derivative_order` and `D2` is an `AbstractDerivativeOperator` of second `derivative_order` or an `AbstractMatrix`.
-Both summation by parts operators should be associated with the same grid.
+Create a solver, where `D1` is an `AbstractDerivativeOperator`
+from [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
+of first `derivative_order` and `D2` is an `AbstractDerivativeOperator`
+of second `derivative_order` or an `AbstractMatrix`. It can also be `nothing`
+if no second derivative is used by the discretization.
+Both summation-by-parts operators should be associated with the same grid.
 """
 function Solver(D1::AbstractDerivativeOperator{RealT},
-                D2::Union{AbstractDerivativeOperator{RealT}, AbstractMatrix{RealT}}) where {
-                                                                                            RealT
-                                                                                            }
+                D2::Union{AbstractDerivativeOperator{RealT}, AbstractMatrix{RealT},
+                          Nothing}) where {RealT}
     Solver{RealT, typeof(D1), typeof(D2)}(D1, D2)
 end
 
@@ -102,8 +106,8 @@ end
 end
 
 function allocate_coefficients(mesh::Mesh1D, equations, solver::AbstractSolver)
-    return ArrayPartition([zeros(real(solver), nnodes(mesh))
-                           for _ in eachvariable(equations)]...)
+    return ArrayPartition(ntuple(_ -> zeros(real(solver), nnodes(mesh)),
+                                 Val(nvariables(equations))))
 end
 
 function compute_coefficients!(q, func, t, mesh::Mesh1D, equations, solver::AbstractSolver)
diff --git a/test/Project.toml b/test/Project.toml
index 9c9acb41..f6bbeb9e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,5 +13,5 @@ ExplicitImports = "1.0.1"
 OrdinaryDiffEq = "6.62"
 Plots = "1.25"
 SparseArrays = "1"
-SummationByPartsOperators = "0.5.52"
+SummationByPartsOperators = "0.5.63"
 Test = "1"
diff --git a/test/runtests.jl b/test/runtests.jl
index 55376424..552a6aee 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,4 +7,5 @@ using Test
     include("test_bbm_bbm_1d.jl")
     include("test_bbm_bbm_variable_bathymetry_1d.jl")
     include("test_svaerd_kalisch_1d.jl")
+    include("test_serre_green_naghdi_1d.jl")
 end
diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl
new file mode 100644
index 00000000..c5bb1e55
--- /dev/null
+++ b/test/test_serre_green_naghdi_1d.jl
@@ -0,0 +1,64 @@
+module TestSerreGreenNaghdi1D
+
+using Test
+using DispersiveShallowWater
+
+include("test_util.jl")
+
+EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d")
+
+@testset "SerreGreenNaghdiEquations1D" begin
+    @trixi_testset "serre_green_naghdi_soliton.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_soliton.jl"),
+                            tspan=(0.0, 0.1),
+                            l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0],
+                            linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0],
+                            cons_error=[0.0, 8.174581012099225e-10, 0.0],
+                            change_waterheight=0.0,
+                            change_entropy_modified=-3.1093350116861984e-11)
+
+        @test_allocations(semi, sol, 550_000)
+    end
+
+    @trixi_testset "serre_green_naghdi_soliton_fourier.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_soliton_fourier.jl"),
+                            tspan=(0.0, 0.1),
+                            l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0],
+                            linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0],
+                            cons_error=[2.842170943040401e-14, 4.627409566637652e-13, 0.0],
+                            change_waterheight=2.842170943040401e-14,
+                            change_entropy_modified=-3.097966327914037e-11)
+
+        @test_allocations(semi, sol, allocs=450_000)
+    end
+
+    @trixi_testset "serre_green_naghdi_soliton_upwind.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_soliton_upwind.jl"),
+                            tspan=(0.0, 0.1),
+                            l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0],
+                            linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0],
+                            cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0],
+                            change_waterheight=4.263256414560601e-14,
+                            change_entropy_modified=-3.1036506698001176e-11)
+
+        @test_allocations(semi, sol, allocs=500_000)
+    end
+
+    @trixi_testset "serre_green_naghdi_soliton_relaxation.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_soliton_relaxation.jl"),
+                            tspan=(0.0, 0.1),
+                            l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0],
+                            linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0],
+                            cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0],
+                            change_waterheight=2.842170943040401e-14,
+                            change_entropy_modified=0.0)
+
+        @test_allocations(semi, sol, allocs=450_000)
+    end
+end
+
+end # module
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 24e3418e..41968101 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -99,7 +99,8 @@ using SparseArrays: sparse, SparseMatrixCSC
     end
 
     @testset "BBMBBMEquations1D" begin
-        equations = @test_nowarn BBMBBMEquations1D(gravity_constant = 9.81, D = 2.0)
+        equations = @test_nowarn @inferred BBMBBMEquations1D(gravity_constant = 9.81,
+                                                             D = 2.0)
         @test_nowarn print(equations)
         @test_nowarn display(equations)
         conversion_functions = [
@@ -117,14 +118,34 @@ using SparseArrays: sparse, SparseMatrixCSC
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
         end
         q = [42.0, 2.0]
-        @test prim2prim(q, equations) == q
+        @test @inferred(prim2prim(q, equations)) == q
         @test isapprox(cons2prim(prim2cons(q, equations), equations), q)
-        @test waterheight_total(q, equations) == 42.0
-        @test waterheight(q, equations) == 44.0
-        @test velocity(q, equations) == 2.0
-        @test momentum(q, equations) == 88.0
-        @test discharge(q, equations) == 88.0
-        @test isapprox(energy_total(q, equations), 8740.42)
+        @test @inferred(waterheight_total(q, equations)) == 42.0
+        @test @inferred(waterheight(q, equations)) == 44.0
+        @test @inferred(velocity(q, equations)) == 2.0
+        @test @inferred(momentum(q, equations)) == 88.0
+        @test @inferred(discharge(q, equations)) == 88.0
+        @test @inferred(still_water_surface(q, equations)) == 0.0
+        @test isapprox(@inferred(energy_total(q, equations)), 8740.42)
+        @test @inferred(energy_total(q, equations)) == @inferred(entropy(q, equations))
+
+        @testset "default implementation of energy_total_modified" begin
+            initial_condition = initial_condition_convergence_test
+            boundary_conditions = boundary_condition_periodic
+            mesh = @inferred Mesh1D(-1.0, 1.0, 10)
+            solver = Solver(mesh, 4)
+            semi = @inferred Semidiscretization(mesh, equations, initial_condition,
+                                                solver; boundary_conditions)
+            q = @inferred DispersiveShallowWater.compute_coefficients(initial_condition,
+                                                                      0.0, semi)
+            _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
+            e_modified = @inferred energy_total_modified(q, equations, cache)
+            e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(q -> energy_total(q,
+                                                                                            equations),
+                                                                          q, semi)
+            @test isapprox(e_modified_total, e_total)
+        end
     end
 
     @testset "BBMBBMVariableEquations1D" begin
@@ -146,14 +167,34 @@ using SparseArrays: sparse, SparseMatrixCSC
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
         end
         q = [42.0, 2.0, 2.0]
-        @test prim2prim(q, equations) == q
-        @test isapprox(cons2prim(prim2cons(q, equations), equations), q)
-        @test waterheight_total(q, equations) == 42.0
-        @test waterheight(q, equations) == 44.0
-        @test velocity(q, equations) == 2.0
-        @test momentum(q, equations) == 88.0
-        @test discharge(q, equations) == 88.0
-        @test isapprox(energy_total(q, equations), 8740.42)
+        @test @inferred(prim2prim(q, equations)) == q
+        @test isapprox(@inferred(cons2prim(prim2cons(q, equations), equations)), q)
+        @test @inferred(waterheight_total(q, equations)) == 42.0
+        @test @inferred(waterheight(q, equations)) == 44.0
+        @test @inferred(velocity(q, equations)) == 2.0
+        @test @inferred(momentum(q, equations)) == 88.0
+        @test @inferred(discharge(q, equations)) == 88.0
+        @test @inferred(still_water_surface(q, equations)) == 0.0
+        @test isapprox(@inferred(energy_total(q, equations)), 8740.42)
+        @test @inferred(energy_total(q, equations)) == @inferred(entropy(q, equations))
+
+        @testset "default implementation of energy_total_modified" begin
+            initial_condition = initial_condition_convergence_test
+            boundary_conditions = boundary_condition_periodic
+            mesh = @inferred Mesh1D(-1.0, 1.0, 10)
+            solver = Solver(mesh, 4)
+            semi = @inferred Semidiscretization(mesh, equations, initial_condition,
+                                                solver; boundary_conditions)
+            q = @inferred DispersiveShallowWater.compute_coefficients(initial_condition,
+                                                                      0.0, semi)
+            _, _, _, cache = @inferred DispersiveShallowWater.mesh_equations_solver_cache(semi)
+            e_modified = @inferred energy_total_modified(q, equations, cache)
+            e_modified_total = @inferred DispersiveShallowWater.integrate(e_modified, semi)
+            e_total = @inferred DispersiveShallowWater.integrate_quantity(q -> energy_total(q,
+                                                                                            equations),
+                                                                          q, semi)
+            @test isapprox(e_modified_total, e_total)
+        end
     end
 
     @testset "SvaerdKalischEquations1D" begin
@@ -180,14 +221,44 @@ using SparseArrays: sparse, SparseMatrixCSC
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
         end
         q = [42.0, 2.0, 2.0]
-        @test prim2prim(q, equations) == q
-        @test isapprox(cons2prim(prim2cons(q, equations), equations), q)
-        @test waterheight_total(q, equations) == 42.0
-        @test waterheight(q, equations) == 44.0
-        @test velocity(q, equations) == 2.0
-        @test momentum(q, equations) == 88.0
-        @test discharge(q, equations) == 88.0
-        @test isapprox(energy_total(q, equations), 8740.42)
+        @test @inferred(prim2prim(q, equations)) == q
+        @test isapprox(@inferred(cons2prim(prim2cons(q, equations), equations)), q)
+        @test @inferred(waterheight_total(q, equations)) == 42.0
+        @test @inferred(waterheight(q, equations)) == 44.0
+        @test @inferred(velocity(q, equations)) == 2.0
+        @test @inferred(momentum(q, equations)) == 88.0
+        @test @inferred(discharge(q, equations)) == 88.0
+        @test @inferred(still_water_surface(q, equations)) == 0.0
+        @test isapprox(@inferred(energy_total(q, equations)), 8740.42)
+    end
+
+    @testset "SerreGreenNaghdiEquations1D" begin
+        equations = @test_nowarn SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+        @test_nowarn print(equations)
+        @test_nowarn display(equations)
+        conversion_functions = [
+            waterheight_total,
+            waterheight,
+            velocity,
+            momentum,
+            discharge,
+            entropy,
+            energy_total,
+            prim2cons,
+            prim2prim,
+        ]
+        for conversion in conversion_functions
+            @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
+        end
+        q = [42.0, 2.0, 0.0]
+        @test @inferred(prim2prim(q, equations)) == q
+        @test isapprox(@inferred(cons2prim(prim2cons(q, equations), equations)), q)
+        @test @inferred(waterheight_total(q, equations)) == 42.0
+        @test @inferred(waterheight(q, equations)) == 42.0
+        @test @inferred(velocity(q, equations)) == 2.0
+        @test @inferred(momentum(q, equations)) == 84.0
+        @test @inferred(discharge(q, equations)) == 84.0
+        @test @inferred(still_water_surface(q, equations)) == 0.0
     end
 
     @testset "AnalysisCallback" begin