diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl
index 3949dbb4..886167f8 100644
--- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl
+++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl
@@ -4,7 +4,7 @@ using DispersiveShallowWater
 ###############################################################################
 # Semidiscretization of the BBM-BBM equations
 
-equations = BBMBBMVariableEquations1D(gravity_constant = 9.81, eta0 = 0.8)
+equations = BBMBBMVariableEquations1D(gravity_constant = 9.81)
 
 initial_condition = initial_condition_dingemans
 boundary_conditions = boundary_condition_periodic
diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl
index ad1e1666..ab640e15 100644
--- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl
+++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl
@@ -6,7 +6,7 @@ using SparseArrays: sparse
 ###############################################################################
 # Semidiscretization of the BBM-BBM equations
 
-equations = BBMBBMVariableEquations1D(gravity_constant = 1.0, eta0 = 2.0)
+equations = BBMBBMVariableEquations1D(gravity_constant = 1.0)
 
 # Setup a truly discontinuous bottom topography function for this academic
 # testcase of well-balancedness. The errors from the analysis callback are
@@ -18,11 +18,11 @@ function initial_condition_discontinuous_well_balancedness(x, t,
     # Set the background values
     eta = equations.eta0
     v = 0.0
-    D = -1.0
+    D = equations.eta0 - 1.0
 
     # Setup a discontinuous bottom topography
     if x >= 0.5 && x <= 0.75
-        D = -1.5 - 0.5 * sinpi(2.0 * x)
+        D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x)
     end
 
     return SVector(eta, v, D)
diff --git a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl
index 759713f0..9f0403fd 100644
--- a/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl
+++ b/examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl
@@ -19,11 +19,11 @@ function initial_condition_discontinuous_well_balancedness(x, t,
     # Set the background values
     eta = equations.eta0
     v = 0.0
-    D = -1.0
+    D = equations.eta0 - 1.0
 
     # Setup a discontinuous bottom topography
     if x >= 0.5 && x <= 0.75
-        D = -1.5 - 0.5 * sinpi(2.0 * x)
+        D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x)
     end
 
     return SVector(eta, v, D)
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 3dc986e6..3e476957 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMEquations1D(gravity, D)
+    BBMBBMEquations1D(gravity, D, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equations are given by
 ```math
@@ -9,8 +9,8 @@ BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equation
 \end{aligned}
 ```
 The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the constant bottom topography (bathymetry) ``b = -D``. The water height above the bathymetry is therefore given by
-``h = \eta + D``.
+The gravitational constant is denoted by `g` and the constant bottom topography (bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry is therefore given by
+``h = \eta - \eta_0 + D``. The BBM-BBM equations are only implemented for ``\eta_0 = 0``.
 
 One reference for the BBM-BBM system can be found in
 - Jerry L. Bona, Min Chen (1998)
@@ -21,10 +21,12 @@ One reference for the BBM-BBM system can be found in
 struct BBMBBMEquations1D{RealT <: Real} <: AbstractBBMBBMEquations{1, 2}
     gravity::RealT # gravitational constant
     D::RealT       # constant bathymetry
+    eta0::RealT    # constant still-water surface
 end
 
-function BBMBBMEquations1D(; gravity_constant, D = 1.0)
-    BBMBBMEquations1D(gravity_constant, D)
+function BBMBBMEquations1D(; gravity_constant, D = 1.0, eta0 = 0.0)
+    eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM-BBM equations"
+    BBMBBMEquations1D(gravity_constant, D, eta0)
 end
 
 varnames(::typeof(prim2prim), ::BBMBBMEquations1D) = ("η", "v")
@@ -94,11 +96,12 @@ function create_cache(mesh,
                       initial_condition,
                       RealT,
                       uEltype)
+    D = equations.D
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
+        invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
     elseif solver.D1 isa PeriodicUpwindOperators
-        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
+        invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
     else
         @error "unknown type of first-derivative operator"
     end
@@ -122,15 +125,15 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond
     deta = view(dq, 1, :)
     dv = view(dq, 2, :)
 
+    D = equations.D
     # energy and mass conservative semidiscretization
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (equations.D * v + eta .* v)
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (D * v + eta .* v)
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1 *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     elseif solver.D1 isa PeriodicUpwindOperators
-        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.central *
-                                                  (equations.D * v + eta .* v)
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.central * (D * v + eta .* v)
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.central *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     else
@@ -148,7 +151,7 @@ end
 @inline function prim2cons(q, equations::BBMBBMEquations1D)
     eta, v = q
 
-    h = eta + equations.D
+    h = eta - equations.eta0 + equations.D
     hv = h * v
     return SVector(h, hv)
 end
@@ -156,7 +159,7 @@ end
 @inline function cons2prim(u, equations::BBMBBMEquations1D)
     h, hv = u
 
-    eta = h - equations.D
+    eta = h + equations.eta0 - equations.D
     v = hv / h
     return SVector(eta, v)
 end
@@ -170,7 +173,7 @@ end
 end
 
 @inline function bathymetry(q, equations::BBMBBMEquations1D)
-    return -equations.D
+    return equations.eta0 - equations.D
 end
 
 @inline function waterheight(q, equations::BBMBBMEquations1D)
@@ -179,7 +182,8 @@ end
 
 @inline function energy_total(q, equations::BBMBBMEquations1D)
     eta, v = q
-    e = 0.5 * (equations.gravity * eta^2 + (equations.D + eta) * v^2)
+    D = still_waterdepth(q, equations)
+    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
     return e
 end
 
diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
index bc8525ef..5019c4bc 100644
--- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl
+++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl
@@ -9,8 +9,8 @@ BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension with spatiall
 \end{aligned}
 ```
 The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = -D``. The water height above the bathymetry is therefore given by
-``h = \eta + D``.
+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 BBM-BBM equations are only implemented for ``\eta_0 = 0``.
 
 One reference for the BBM-BBM system with spatially varying bathymetry can be found in
 - Samer Israwi, Henrik Kalisch, Theodoros Katsaounis, Dimitrios Mitsotakis (2022)
@@ -20,10 +20,11 @@ One reference for the BBM-BBM system with spatially varying bathymetry can be fo
 """
 struct BBMBBMVariableEquations1D{RealT <: Real} <: AbstractBBMBBMEquations{1, 3}
     gravity::RealT # gravitational constant
-    eta0::RealT    # constant "lake-at-rest" total water height
+    eta0::RealT    # constant still-water surface
 end
 
-function BBMBBMVariableEquations1D(; gravity_constant, eta0 = 1.0)
+function BBMBBMVariableEquations1D(; gravity_constant, eta0 = 0.0)
+    eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM-BBM equations"
     BBMBBMVariableEquations1D(gravity_constant, eta0)
 end
 
@@ -107,6 +108,11 @@ The initial condition that uses the dispersion relation of the Euler equations
 to approximate waves generated by a wave maker as it is done by experiments of
 Dingemans. The topography is a trapezoidal.
 
+!!! warning "Translation of water height"
+    The initial condition for the water height is translated to be around 0, which is
+    needed for the simulation because the `BBMBBMVariableEquations1D` are only implemented
+    for ``\eta_0 = 0`.
+
 References:
 - Magnus Svärd, Henrik Kalisch (2023)
   A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
@@ -116,16 +122,16 @@ References:
   [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
 """
 function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh)
-    eta0 = 0.8
+    h0 = 0.8
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
-    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl
+    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl
     if x < -30.5 * pi / k || x > -8.5 * pi / k
         h = 0.0
     else
         h = A * cos(k * x)
     end
-    v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0
+    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
     if 11.01 <= x && x < 23.04
         b = 0.6 * (x - 11.01) / (23.04 - 11.01)
     elseif 23.04 <= x && x < 27.04
@@ -135,8 +141,10 @@ function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D,
     else
         b = 0.0
     end
-    eta = h + eta0
-    D = -b
+    # Here, we compute eta - h0!! To obtain the original eta, h0 = 0.8 needs to be added again!
+    # This is because the BBM-BBM equations are only implemented for eta0 = 0
+    eta = h
+    D = h0 - b
     return SVector(eta, v, D)
 end
 
@@ -150,7 +158,7 @@ function create_cache(mesh,
     D = Array{RealT}(undef, nnodes(mesh))
     x = grid(solver)
     for i in eachnode(solver)
-        D[i] = initial_condition(x[i], 0.0, equations, mesh)[3]
+        D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)
     end
     K = Diagonal(D .^ 2)
     if solver.D1 isa PeriodicDerivativeOperator ||
@@ -164,7 +172,7 @@ function create_cache(mesh,
         @error "unknown type of first-derivative operator"
     end
     tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
-    return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1)
+    return (invImDKD = invImDKD, invImD2K = invImD2K, D = D, tmp1 = tmp1)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -175,14 +183,13 @@ end
 function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack invImDKD, invImD2K = cache
+    @unpack invImDKD, invImD2K, D = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
 
     eta = view(q, 1, :)
     v = view(q, 2, :)
-    D = view(q, 3, :)
     deta = view(dq, 1, :)
     dv = view(dq, 2, :)
     dD = view(dq, 3, :)
@@ -212,9 +219,9 @@ end
 @inline function prim2cons(q, equations::BBMBBMVariableEquations1D)
     eta, v, D = q
 
-    h = eta + D
+    b = bathymetry(q, equations)
+    h = eta - b
     hv = h * v
-    b = -D
     return SVector(h, hv, b)
 end
 
@@ -223,7 +230,7 @@ end
 
     eta = h + b
     v = hv / h
-    D = -b
+    D = equations.eta0 - b
     return SVector(eta, v, D)
 end
 
@@ -236,7 +243,8 @@ end
 end
 
 @inline function bathymetry(q, equations::BBMBBMVariableEquations1D)
-    return -q[3]
+    D = q[3]
+    return equations.eta0 - D
 end
 
 @inline function waterheight(q, equations::BBMBBMVariableEquations1D)
@@ -245,7 +253,7 @@ end
 
 @inline function energy_total(q, equations::BBMBBMVariableEquations1D)
     eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta) * v^2)
+    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
     return e
 end
 
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index ba3d02b3..ba828f92 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -133,6 +133,12 @@ See [`momentum`](@ref).
 
 varnames(::typeof(discharge), equations) = ("P",)
 
+@inline function still_waterdepth(q, equations::AbstractEquations)
+    b = bathymetry(q, equations)
+    D = equations.eta0 - b
+    return D
+end
+
 """
     entropy(q, equations)
 
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index f537ab13..00e33777 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -8,7 +8,7 @@ Dispersive system by Svärd and Kalisch in one spatial dimension with spatially
   (hv)_t + (hv^2)_x + gh(h + b)_x &= (\hat\alpha v(\hat\alpha(h + b)_x)_x)_x + (\hat\beta v_x)_{xt} + \frac{1}{2}(\hat\gamma v_x)_{xx} + \frac{1}{2}(\hat\gamma v_{xx})_x,
 \end{aligned}
 ```
-where ``\hat\alpha^2 = \alpha\sqrt{gd}d^2``, ``\hat\beta = \beta d^3``, ``\hat\gamma = \gamma\sqrt{gd}d^3``. The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``d = \eta_0 - b`` is the still-water depth and `eta0` is the still-water surface (lake-at-rest).
+where ``\hat\alpha^2 = \alpha\sqrt{gD}D^2``, ``\hat\beta = \beta D^3``, ``\hat\gamma = \gamma\sqrt{gD}D^3``. The coefficients ``\alpha``, ``\beta`` and ``\gamma`` are provided in dimensionless form and ``D = \eta_0 - b`` is the still-water depth and `eta0` is the still-water surface (lake-at-rest).
 The equations can be rewritten in primitive variables as
 ```math
 \begin{aligned}
@@ -17,8 +17,8 @@ The equations can be rewritten in primitive variables as
 \end{aligned}
 ```
 The unknown quantities of the Svärd-Kalisch equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = -D``. The water height above the bathymetry is therefore given by
-``h = \eta + D``.
+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 equations by Svärd and Kalisch are presented and analyzed in
 - Magnus Svärd, Henrik Kalisch (2023)
@@ -28,7 +28,7 @@ The equations by Svärd and Kalisch are presented and analyzed in
 """
 struct SvaerdKalischEquations1D{RealT <: Real} <: AbstractSvaerdKalischEquations{1, 3}
     gravity::RealT # gravitational constant
-    eta0::RealT    # constant "lake-at-rest" total water height
+    eta0::RealT    # constant still-water surface
     alpha::RealT   # coefficient
     beta::RealT    # coefficient
     gamma::RealT   # coefficient
@@ -36,7 +36,7 @@ end
 
 const SvärdKalischEquations1D = SvaerdKalischEquations1D
 
-function SvaerdKalischEquations1D(; gravity_constant, eta0 = 1.0, alpha = 0.0,
+function SvaerdKalischEquations1D(; gravity_constant, eta0 = 0.0, alpha = 0.0,
                                   beta = 0.2308939393939394, gamma = 0.04034343434343434)
     SvaerdKalischEquations1D(gravity_constant, eta0, alpha, beta, gamma)
 end
@@ -60,16 +60,16 @@ References:
   [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
 """
 function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D, mesh)
-    eta0 = 0.8
+    h0 = 0.8
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
-    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl
+    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl
     if x < -30.5 * pi / k || x > -8.5 * pi / k
         h = 0.0
     else
         h = A * cos(k * x)
     end
-    v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0
+    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
     if 11.01 <= x && x < 23.04
         b = 0.6 * (x - 11.01) / (23.04 - 11.01)
     elseif 23.04 <= x && x < 27.04
@@ -79,8 +79,8 @@ function initial_condition_dingemans(x, t, equations::SvaerdKalischEquations1D,
     else
         b = 0.0
     end
-    eta = h + eta0
-    D = -b
+    eta = h + h0
+    D = equations.eta0 - b
     return SVector(eta, v, D)
 end
 
@@ -170,14 +170,13 @@ function create_cache(mesh,
     D = Array{RealT}(undef, nnodes(mesh))
     x = grid(solver)
     for i in eachnode(solver)
-        D[i] = initial_condition(x[i], 0.0, equations, mesh)[3]
+        D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)
     end
-    d = equations.eta0 .+ D
     h = Array{RealT}(undef, nnodes(mesh))
     hv = similar(h)
-    alpha_hat = sqrt.(equations.alpha * sqrt.(equations.gravity * d) .* d .^ 2)
-    beta_hat = equations.beta * d .^ 3
-    gamma_hat = equations.gamma * sqrt.(equations.gravity * d) .* d .^ 3
+    alpha_hat = sqrt.(equations.alpha * sqrt.(equations.gravity * D) .* D .^ 2)
+    beta_hat = equations.beta * D .^ 3
+    gamma_hat = equations.gamma * sqrt.(equations.gravity * D) .* D .^ 3
     tmp1 = similar(h)
     tmp2 = similar(h)
     hmD1betaD1 = Array{RealT}(undef, nnodes(mesh), nnodes(mesh))
@@ -192,7 +191,7 @@ function create_cache(mesh,
     else
         @error "unknown type of first-derivative operator"
     end
-    return (hmD1betaD1 = hmD1betaD1, D1betaD1 = D1betaD1, d = d, h = h, hv = hv,
+    return (hmD1betaD1 = hmD1betaD1, D1betaD1 = D1betaD1, D = D, h = h, hv = hv,
             alpha_hat = alpha_hat, beta_hat = beta_hat, gamma_hat = gamma_hat,
             tmp1 = tmp1, tmp2 = tmp2, D1_central = D1_central, D1 = solver.D1)
 end
@@ -201,20 +200,19 @@ end
 function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack hmD1betaD1, D1betaD1, d, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache
+    @unpack hmD1betaD1, D1betaD1, D, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
 
     eta = view(q, 1, :)
     v = view(q, 2, :)
-    D = view(q, 3, :)
     deta = view(dq, 1, :)
     dv = view(dq, 2, :)
     dD = view(dq, 3, :)
     fill!(dD, zero(eltype(dD)))
 
     @timeit timer() "deta hyperbolic" begin
-        h = eta .+ D
+        h = eta .+ D .- equations.eta0
         hv = h .* v
 
         if solver.D1 isa PeriodicDerivativeOperator ||
@@ -268,9 +266,9 @@ end
 @inline function prim2cons(q, equations::SvaerdKalischEquations1D)
     eta, v, D = q
 
-    h = eta + D
+    b = bathymetry(q, equations)
+    h = eta - b
     hv = h * v
-    b = -D
     return SVector(h, hv, b)
 end
 
@@ -279,7 +277,7 @@ end
 
     eta = h + b
     v = hv / h
-    D = -b
+    D = equations.eta0 - b
     return SVector(eta, v, D)
 end
 
@@ -292,7 +290,8 @@ end
 end
 
 @inline function bathymetry(q, equations::SvaerdKalischEquations1D)
-    return -q[3]
+    D = q[3]
+    return equations.eta0 - D
 end
 
 @inline function waterheight(q, equations::SvaerdKalischEquations1D)
@@ -301,7 +300,7 @@ end
 
 @inline function energy_total(q, equations::SvaerdKalischEquations1D)
     eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta) * v^2)
+    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
     return e
 end
 
@@ -327,7 +326,7 @@ number of nodes as length of the second dimension.
     eta = view(q, 1, :)
     v = view(q, 2, :)
     D = view(q, 3, :)
-    beta_hat = equations.beta * (eta .+ D) .^ 3
+    beta_hat = equations.beta * D .^ 3
     if cache.D1 isa PeriodicDerivativeOperator ||
        cache.D1 isa UniformPeriodicCoupledOperator
         tmp = 0.5 * beta_hat .* ((cache.D1 * v) .^ 2)
diff --git a/test/test_bbm_bbm_variable_bathymetry_1d.jl b/test/test_bbm_bbm_variable_bathymetry_1d.jl
index 7910ab09..03f2abf2 100644
--- a/test/test_bbm_bbm_variable_bathymetry_1d.jl
+++ b/test/test_bbm_bbm_variable_bathymetry_1d.jl
@@ -90,12 +90,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_variable_bathymetry_1d")
                                      "bbm_bbm_variable_bathymetry_1d_dingemans.jl"),
                             tspan=(0.0, 1.0),
                             N=512,
-                            l2=[0.2322073114427607 0.7753458687737584 0.0],
-                            linf=[0.037222719015511885 0.124336213226626 0.0],
+                            l2=[0.2229215734522602 0.750492441160796 0.0],
+                            linf=[0.03584030574058169 0.1202292994661 0.0],
                             cons_error=[1.4210854715202004e-13 3.1478183774857893e-15 0.0],
                             change_waterheight=-1.4210854715202004e-13,
                             change_velocity=-3.1478183774857893e-15,
-                            change_entropy=-1.442231223336421e-9)
+                            change_entropy=3.417533493699221e-7)
     end
 end
 
diff --git a/test/test_svaerd_kalisch_1d.jl b/test/test_svaerd_kalisch_1d.jl
index 979d5a05..f120f3f2 100644
--- a/test/test_svaerd_kalisch_1d.jl
+++ b/test/test_svaerd_kalisch_1d.jl
@@ -30,7 +30,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             cons_error=[3.979039320256561e-13 4.937137540373564e-5 0.0],
                             change_waterheight=-3.979039320256561e-13,
                             change_entropy=-0.00024362648639453255,
-                            change_entropy_modified=3.240868750253867e-6)
+                            change_entropy_modified=-6.311893230304122e-9)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_cg" begin
@@ -42,7 +42,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             cons_error=[1.4210854715202004e-13 4.953054817174909e-5 0.0],
                             change_waterheight=-1.4210854715202004e-13,
                             change_entropy=-0.0002425303440531934,
-                            change_entropy_modified=2.9729466177741415e-6)
+                            change_entropy_modified=-2.6815314413397573e-9)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_upwind" begin
@@ -54,7 +54,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                             cons_error=[1.1368683772161603e-13 4.871598417571836e-5 0.0],
                             change_waterheight=-1.1368683772161603e-13,
                             change_entropy=-0.00023645232727176335,
-                            change_entropy_modified=2.0682244894487667e-6)
+                            change_entropy_modified=-6.654090611846186e-9)
     end
 
     @trixi_testset "svaerd_kalisch_1d_dingemans_relaxation" begin
@@ -62,11 +62,11 @@ EXAMPLES_DIR = joinpath(examples_dir(), "svaerd_kalisch_1d")
                                      "svaerd_kalisch_1d_dingemans_relaxation.jl"),
                             tspan=(0.0, 1.0),
                             N=512,
-                            l2=[0.2279601359992612 0.7519295302772953 0.0],
-                            linf=[0.03670816991272552 0.12141112705223758 0.0],
-                            cons_error=[3.979039320256561e-13 4.937017087502672e-5 0.0],
+                            l2=[0.22796107242561717 0.7519327155905444 0.0],
+                            linf=[0.03670834831604197 0.12141172368792873 0.0],
+                            cons_error=[3.979039320256561e-13 4.937137655207271e-5 0.0],
                             change_waterheight=-3.979039320256561e-13,
-                            change_entropy=-0.0002466715038735856,
+                            change_entropy=-0.00024362054875837202,
                             change_entropy_modified=0.0)
     end
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 17130564..f641cdf8 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -108,7 +108,7 @@ using SparseArrays: sparse, SparseMatrixCSC
             entropy,
             energy_total,
             prim2cons,
-            prim2prim,
+            prim2prim
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -137,7 +137,7 @@ using SparseArrays: sparse, SparseMatrixCSC
             entropy,
             energy_total,
             prim2cons,
-            prim2prim,
+            prim2prim
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple
@@ -171,7 +171,7 @@ using SparseArrays: sparse, SparseMatrixCSC
             prim2cons,
             prim2prim,
             energy_total_modified,
-            entropy_modified,
+            entropy_modified
         ]
         for conversion in conversion_functions
             @test DispersiveShallowWater.varnames(conversion, equations) isa Tuple