diff --git a/docs/src/basic_example.md b/docs/src/basic_example.md
index 87515b7c..16d8cd51 100644
--- a/docs/src/basic_example.md
+++ b/docs/src/basic_example.md
@@ -131,7 +131,7 @@ By default, this will plot the bathymetry, but not the initial (analytical) solu
 
 You can adjust this by passing the boolean values `plot_bathymetry` (if `true`, always plot bathymetry in the first subplot) and `plot_initial`. Note that `plot_initial = true` will evaluate and plot the initial condition function at the same time `t` as the numerical solution being displayed (the final time by default). This means if your initial condition function represents an analytical solution, setting `plot_initial = true` will plot the analytical solution at that specific time for comparison.
 
-Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step. Here `conversion = waterheight_total` makes it so that we only look at the waterheight ``\eta`` and not also the velocity ``v``. More on tutorials for plotting can be found in the chapter [Plotting Simulation Results](@ref plotting).
+Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step. Here `conversion = waterheight_total` makes it so that we only look at the total water height ``\eta`` and not also the velocity ``v``. More on tutorials for plotting can be found in the chapter [Plotting Simulation Results](@ref plotting).
 
 ```@example overview
 anim = @animate for step in 1:length(sol.u)
diff --git a/docs/src/callbacks.md b/docs/src/callbacks.md
index 3070fa64..3ed463ca 100644
--- a/docs/src/callbacks.md
+++ b/docs/src/callbacks.md
@@ -77,13 +77,17 @@ nothing # hide
 
 The analysis callback computes ``L^2`` and ``L^\infty`` errors by comparing the numerical solution to the initial condition at time ``t`` (which can be the analytical solution, if available). Additional error types can be specified using the `extra_analysis_errors` parameter, and physical quantities can be monitored using `extra_analysis_integrals`.
 
-The conservation error measures the temporal change of conserved quantities. For the BBM-BBM equations, important conserved quantities include the total water mass (integral of water height `h`), the total momentum (integral of `v` for flat bathymetry), and the total [`entropy`](@ref). The specific form of the entropy varies between different equation systems. For the BBM-BBM equations, the integral of the entropy is:
+The conservation error measures the temporal change of conserved quantities. For the BBM-BBM equations, important conserved quantities include the total water mass (integral of water height ``h`` given
+by [`waterheight`](@ref)), the total momentum (integral of ``v`` for flat bathymetry given by [`velocity`](@ref)), the integral of the total water height ``\eta`` given by [`waterheight_total`](@ref), and
+the total entropy (integral of ``U`` given by [`entropy`](@ref)). The specific form of the [`entropy`](@ref) varies between different equation systems. For the BBM-BBM equations, the integral of the entropy is:
 
 ```math
-\mathcal E(t; \eta, v) = \frac{1}{2}\int_\Omega g\eta^2 + (\eta + D)v^2 \, dx
+\mathcal E(t; \eta, v) = \int_\Omega U(\eta, v) \, dx = \frac{1}{2}\int_\Omega g\eta^2 + (\eta + D)v^2 \, dx
 ```
 
-where ``\eta`` is the total water height and ``D`` is the still-water depth.
+where ``\eta`` is the total water height and ``D`` is the still-water depth. Any scalar quantity taking a vector `q` of primitive variables and the `equations` as arguments and returning a scalar value can be used
+as an additional integral to monitor during the simulation by passing it to the `extra_analysis_integrals` parameter. Other quantities that can be monitored include [`momentum`](@ref), [`entropy_modified`](@ref),
+and [`hamiltonian`](@ref).
 
 ```@example callback
 tspan = (0.0, 20.0)
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 4d59848d..7f044c03 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -336,6 +336,7 @@ function analyze!(semi::Semidiscretization, quantity, q, t)
     integrate_quantity!(semi.cache.tmp1, quantity, q, semi)
 end
 
+pretty_form_utf(::typeof(waterheight)) = "∫h"
 pretty_form_utf(::typeof(waterheight_total)) = "∫η"
 pretty_form_utf(::typeof(velocity)) = "∫v"
 pretty_form_utf(::typeof(momentum)) = "∫P"
diff --git a/src/equations/bbm_1d.jl b/src/equations/bbm_1d.jl
index d9282ea1..6fd4b672 100644
--- a/src/equations/bbm_1d.jl
+++ b/src/equations/bbm_1d.jl
@@ -22,12 +22,12 @@ for `split_form = true` and in Linders, Ranocha, and Birken (2023) for `split_fo
 If `split_form` is `true`, a split form in the semidiscretization is used, which conserves
 - the total water mass (integral of ``h``) as a linear invariant
 - a quadratic invariant (integral of ``1/2\eta(\eta - 1/6D^2\eta_{xx})`` or for periodic boundary conditions
-  equivalently ``1/2(\eta^2 + 1/6D^2\eta_x^2)``), which is called here [`energy_total_modified`](@ref)
+  equivalently ``\hat U = 1/2(\eta^2 + 1/6D^2\eta_x^2)``), which is called here [`energy_total_modified`](@ref)
   (and [`entropy_modified`](@ref)) because it contains derivatives of the solution
 
 for periodic boundary conditions. If `split_form` is `false` the semidiscretization conserves
 - the total water mass (integral of ``h``) as a linear invariant
-- the Hamiltonian (integral of ``1/4\sqrt{g/D}\eta^3 + 1/2\sqrt{gD}\eta^2``) (see [`hamiltonian`](@ref))
+- the Hamiltonian (integral of ``H = 1/4\sqrt{g/D}\eta^3 + 1/2\sqrt{gD}\eta^2``) (see [`hamiltonian`](@ref))
 
 for periodic boundary conditions.
 
@@ -195,15 +195,15 @@ function rhs!(dq, q, t, mesh, equations::BBMEquation1D, initial_condition,
 end
 
 """
-    energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
+    energy_total_modified!(e_mod, q_global, equations::BBMEquation1D, cache)
 
-Return the modified total energy `e` of the primitive variables `q_global` for the
+Return the modified total energy ``\\hat e`` of the primitive variables `q_global` for the
 [`BBMEquation1D`](@ref). The `energy_total_modified`
 is a conserved quantity (for periodic boundary conditions).
 
 It is given by
 ```math
-\\frac{1}{2} \\eta(\\eta - \\frac{1}{6}D^2\\eta_{xx}).
+\\hat e(\\eta) = \\frac{1}{2} \\eta(\\eta - \\frac{1}{6}D^2\\eta_{xx}).
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
@@ -211,7 +211,7 @@ It is given by
 
 See also [`energy_total_modified`](@ref).
 """
-function energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
+function energy_total_modified!(e_mod, q_global, equations::BBMEquation1D, cache)
     eta, = q_global.x
     D = equations.D
 
@@ -223,17 +223,17 @@ function energy_total_modified!(e, q_global, equations::BBMEquation1D, cache)
         mul!(eta_xx, D2, eta)
     end
 
-    @.. e = 0.5 * eta * (eta - 1 / 6 * D^2 * eta_xx)
-    return e
+    @.. e_mod = 0.5 * eta * (eta - 1 / 6 * D^2 * eta_xx)
+    return e_mod
 end
 
 """
     hamiltonian!(H, q_global, equations::BBMEquation1D, cache)
 
-Return the Hamiltonian `H` of the primitive variables `q_global` for the
+Return the Hamiltonian ``H`` of the primitive variables `q_global` for the
 [`BBMEquation1D`](@ref). The Hamiltonian is given by
 ```math
-\\frac{1}{4}\\sqrt{\\frac{g}{D}}\\eta^3 + \\frac{1}{2}\\sqrt{gD}\\eta^2.
+H(\\eta) = \\frac{1}{4}\\sqrt{\\frac{g}{D}}\\eta^3 + \\frac{1}{2}\\sqrt{gD}\\eta^2.
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl
index 95ca2e60..f7a70c2f 100644
--- a/src/equations/bbm_bbm_1d.jl
+++ b/src/equations/bbm_bbm_1d.jl
@@ -32,12 +32,12 @@ Ranocha et al. (2020) and generalized for a variable bathymetry in
 Lampert and Ranocha (2025). It conserves
 - the total water mass (integral of ``h``) as a linear invariant
 - the total velocity (integral of ``v``) as a linear invariant for flat bathymetry
-- the total energy
+- the total entropy/energy (integral of ``U = 1/2 hv^2 + 1/2 g\eta^2``)
 
 for periodic boundary conditions (see Lampert, Ranocha). For reflecting boundary conditions,
 the semidiscretization conserves
 - the total water (integral of ``h``) as a linear invariant
-- the total energy.
+- the total entropy/energy (integral of ``U = 1/2 hv^2 + 1/2 g\eta^2``)
 
 Additionally, it is well-balanced for the lake-at-rest stationary solution, see Lampert and Ranocha (2025).
 
@@ -412,7 +412,7 @@ end
 # Discretization that conserves
 # - the total water (integral of ``h``) as a linear invariant
 # - the total momentum (integral of ``v``) as a linear invariant for flat bathymetry
-# - the total energy
+# - the total entropy/energy (integral of ``U = 1/2 hv^2 + 1/2 g\eta^2``)
 # for periodic boundary conditions, see
 # - Joshua Lampert, Hendrik Ranocha (2025)
 #   Structure-preserving numerical methods for two nonlinear systems of dispersive wave equations
@@ -482,7 +482,7 @@ end
 
 # Discretization that conserves
 # - the total water (integral of ``h``) as a linear invariant
-# - the total energy
+# - the total entropy/energy (integral of ``U = 1/2 hv^2 + 1/2 g\eta^2``)
 # for reflecting boundary conditions, see
 # - Joshua Lampert, Hendrik Ranocha (2025)
 #   Structure-preserving numerical methods for two nonlinear systems of dispersive wave equations
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index c9a1bb52..82b85e58 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -142,7 +142,7 @@ end
 """
     waterheight_total(q, equations)
 
-Return the total waterheight of the primitive variables `q` for a given set of
+Return the total water height ``\\eta`` of the primitive variables `q` for a given set of
 `equations`, i.e., the [`waterheight`](@ref) ``h`` plus the
 [`bathymetry`](@ref) ``b``.
 
@@ -158,8 +158,8 @@ 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 ``h`` above the bathymetry ``b``.
+Return the water height of the primitive variables `q` for a given set of
+`equations`, i.e., the water height ``h`` above the bathymetry ``b``.
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
 of the correct length `nvariables(equations)`.
@@ -175,7 +175,7 @@ varnames(::typeof(waterheight), equations) = ("h",)
 """
     velocity(q, equations)
 
-Return the velocity of the primitive variables `q` for a given set of
+Return the velocity ``v`` 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
@@ -190,7 +190,7 @@ varnames(::typeof(velocity), equations) = ("v",)
 """
     momentum(q, equations)
 
-Return the momentum/discharge of the primitive variables `q` for a given set of
+Return the momentum/discharge ``hv`` 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
@@ -230,7 +230,7 @@ still_waterdepth(q, equations::AbstractEquations{1, 1}) = equations.D
 """
     bathymetry(q, equations)
 
-Return the bathymetry of the primitive variables `q` for a given set of
+Return the bathymetry ``b`` 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
@@ -268,12 +268,12 @@ end
 """
     energy_total(q, equations)
 
-Return the total energy of the primitive variables `q` for a given set of
+Return the total energy ``e`` of the primitive variables `q` for a given set of
 `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
+e(\\eta, v) = \\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),
@@ -294,7 +294,7 @@ varnames(::typeof(energy_total), equations) = ("e_total",)
 """
     entropy(q, equations)
 
-Return the entropy of the primitive variables `q` for a given set of
+Return the entropy ``U`` of the primitive variables `q` for a given set of
 `equations`. For all [`AbstractShallowWaterEquations`](@ref), the `entropy`
 is just the [`energy_total`](@ref).
 
@@ -311,7 +311,7 @@ varnames(::typeof(entropy), equations) = ("U",)
 """
     energy_total_modified(q_global, equations, cache)
 
-Return the modified total energy of the primitive variables `q_global` for the
+Return the modified total energy ``\\hat e`` 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
@@ -331,25 +331,25 @@ Internally, this function allocates a vector for the output and
 calls [`DispersiveShallowWater.energy_total_modified!`](@ref).
 """
 function energy_total_modified(q_global, equations::AbstractEquations, cache)
-    e = similar(q_global.x[begin])
-    return energy_total_modified!(e, q_global, equations, cache)
+    e_mod = similar(q_global.x[begin])
+    return energy_total_modified!(e_mod, q_global, equations, cache)
 end
 
 """
-    energy_total_modified!(e, q_global, equations, cache)
+    energy_total_modified!(e_mod, q_global, equations, cache)
 
 In-place version of [`energy_total_modified`](@ref).
 """
-function energy_total_modified!(e, q_global, equations::AbstractEquations,
+function energy_total_modified!(e_mod, q_global, equations::AbstractEquations,
                                 cache)
     # `q_global` is an `ArrayPartition` of the primitive variables at all nodes
     @assert nvariables(equations) == length(q_global.x)
 
     for i in eachindex(q_global.x[begin])
-        e[i] = energy_total(get_node_vars(q_global, equations, i), equations)
+        e_mod[i] = energy_total(get_node_vars(q_global, equations, i), equations)
     end
 
-    return e
+    return e_mod
 end
 
 varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
@@ -360,19 +360,19 @@ varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
 Alias for [`energy_total_modified`](@ref).
 """
 @inline function entropy_modified(q_global, equations::AbstractEquations, cache)
-    e = similar(q_global.x[begin])
-    return entropy_modified!(e, q_global, equations, cache)
+    U_mod = similar(q_global.x[begin])
+    return entropy_modified!(U_mod, q_global, equations, cache)
 end
 
 """
-    entropy_modified!(e, q_global, equations, cache)
+    entropy_modified!(U_mod, q_global, equations, cache)
 
 In-place version of [`entropy_modified`](@ref).
 """
-@inline entropy_modified!(e, q_global, equations, cache) = energy_total_modified!(e,
-                                                                                  q_global,
-                                                                                  equations,
-                                                                                  cache)
+@inline entropy_modified!(U_mod, q_global, equations, cache) = energy_total_modified!(U_mod,
+                                                                                      q_global,
+                                                                                      equations,
+                                                                                      cache)
 
 varnames(::typeof(entropy_modified), equations) = ("U_modified",)
 
@@ -380,7 +380,7 @@ varnames(::typeof(entropy_modified), equations) = ("U_modified",)
 """
     hamiltonian(q_global, equations, cache)
 
-Return the Hamiltonian of the primitive variables `q_global` for the
+Return the Hamiltonian ``H`` of the primitive variables `q_global` for the
 `equations`. The Hamiltonian is a conserved quantity and may contain
 derivatives of the solution.
 
diff --git a/src/equations/hyperbolic_serre_green_naghdi_1d.jl b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
index 8705b837..57ddc4ac 100644
--- a/src/equations/hyperbolic_serre_green_naghdi_1d.jl
+++ b/src/equations/hyperbolic_serre_green_naghdi_1d.jl
@@ -74,7 +74,8 @@ References for the hyperbolized Serre-Green-Naghdi system can be found in
 
 The semidiscretization implemented here conserves
 - the total water mass (integral of ``h``) as a linear invariant
-- the total modified energy
+- the total modified entropy/energy (integral of ``\hat U``/``\hat e``), which is called here
+  [`entropy_modified`](@ref) and [`energy_total_modified`](@ref)
 
 for periodic boundary conditions (see Ranocha and Ricchiuto (2025)).
 Additionally, it is well-balanced for the lake-at-rest stationary solution, see
@@ -357,7 +358,7 @@ end
 
 # Discretization that conserves
 # - the total water mass (integral of ``h``) as a linear invariant
-# - the total modified energy
+# - the total modified entropy/energy (integral of ``\hat U``/``\hat e``) as a nonlinear invariant
 # for periodic boundary conditions, see
 # - Hendrik Ranocha and Mario Ricchiuto (2025)
 #   Structure-Preserving Approximations of the Serre-Green-Naghdi
@@ -520,9 +521,9 @@ end
 
 # The entropy/energy takes the whole `q` for every point in space
 """
-    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::HyperbolicSerreGreenNaghdiEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e_mod, q_global, equations::HyperbolicSerreGreenNaghdiEquations1D, cache)
 
-Return the modified total energy `e` of the primitive variables `q_global` for the
+Return the modified total energy ``\\hat e`` of the primitive variables `q_global` for the
 [`HyperbolicSerreGreenNaghdiEquations1D`](@ref).
 It contains additional terms compared to the usual [`energy_total`](@ref)
 modeling non-hydrostatic contributions. The `energy_total_modified`
@@ -531,7 +532,7 @@ is a conserved quantity (for periodic boundary conditions).
 For a [`bathymetry_mild_slope`](@ref) (and a [`bathymetry_flat`](@ref)),
 the total modified energy is given by
 ```math
-\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 +
+\\hat e(\\eta, v, w) = \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 +
 \\frac{1}{6} h w^2 + \\frac{\\lambda}{6} h (1 - \\eta / h)^2.
 ```
 
@@ -539,7 +540,7 @@ the total modified energy is given by
 
 See also [`energy_total_modified`](@ref).
 """
-function energy_total_modified!(e, q_global,
+function energy_total_modified!(e_mod, q_global,
                                 equations::HyperbolicSerreGreenNaghdiEquations1D,
                                 cache)
     # unpack physical parameters and SBP operator `D1`
@@ -555,8 +556,8 @@ function energy_total_modified!(e, q_global,
     @.. h = eta - b
 
     # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 + λ/6 h (1 - H/h)^2
-    @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * w^2 +
-            lambda / 6 * h * (1 - H / h)^2
+    @.. e_mod = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * w^2 +
+                lambda / 6 * h * (1 - H / h)^2
 
-    return e
+    return e_mod
 end
diff --git a/src/equations/kdv_1d.jl b/src/equations/kdv_1d.jl
index 6b6b0307..301fd9f4 100644
--- a/src/equations/kdv_1d.jl
+++ b/src/equations/kdv_1d.jl
@@ -20,10 +20,10 @@ The KdV equation is first introduced by Joseph Valentin Boussinesq (1877) and re
 
 The semidiscretization implemented here is a modification of the one proposed by Biswas, Ketcheson, Ranocha, and Schütz (2025) for the non-dimensionalized KdV equation ``u_t + u u_x + u_{x x x} = 0.``
 
-The semidiscretization looks the following:
+The semidiscretization is given by:
 ```math
 \begin{aligned}
-  \eta_t+\sqrt{g D} D_1\eta+ 1 / 2 \sqrt{g / D} \eta D_1 \eta +  1 / 2 \sqrt{g / D} D_1 \eta^2 +1 / 6 \sqrt{g D} D^2 D_3\eta &= 0.
+  \eta_t + \sqrt{g D} D_1\eta + 1 / 2 \sqrt{g / D} \eta D_1 \eta +  1 / 2 \sqrt{g / D} D_1 \eta^2 + 1 / 6 \sqrt{g D} D^2 D_3\eta &= 0.
 \end{aligned}
 ```
 where ``D_1`` is a first-derivative operator, ``D_3`` a third-derivative operator, and ``D`` the still-water depth.
@@ -31,7 +31,7 @@ where ``D_1`` is a first-derivative operator, ``D_3`` a third-derivative operato
 It conserves
 - the total water mass (integral of ``\eta``) as a linear invariant
 and if upwind operators (``D_3 = D_{1,+} D_1 D_{1,-}``) or wide-stencil operators (``D_3 = D_1^3``) are used for the third derivative, it also conserves
-- the energy (integral of ``1/2\eta^2``)
+- the entropy/total energy (integral of ``U = 1/2\eta^2``)
 
 for periodic boundary conditions.
 
@@ -279,11 +279,11 @@ end
 """
     energy_total(q, equations::KdVEquation1D)
 
-Return the total energy of the primitive variables `q` for the
+Return the total energy ``e`` of the primitive variables `q` for the
 [`KdVEquation1D`](@ref). For the KdV equation, the total energy consists
 only of the potential energy, given by
 ```math
-\\frac{1}{2} g \\eta^2
+e(\\eta) = \\frac{1}{2} g \\eta^2
 ```
 where ``\\eta`` is the [`waterheight_total`](@ref) and ``g`` is the
 [`gravity`](@ref).
@@ -299,7 +299,7 @@ end
 """
     entropy(q, equations::KdVEquation1D)
 
-Return the entropy of the primitive variables `q` for the [`KdVEquation1D`](@ref).
+Return the entropy ``U`` of the primitive variables `q` for the [`KdVEquation1D`](@ref).
 For the KdV equation, the `entropy` is the same as the [`energy_total`](@ref).
 
 `q` is a vector of the primitive variables at a single node, i.e., a vector
diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl
index 93695ca0..22bcdd46 100644
--- a/src/equations/serre_green_naghdi_1d.jl
+++ b/src/equations/serre_green_naghdi_1d.jl
@@ -59,7 +59,8 @@ References for the Serre-Green-Naghdi system can be found in
 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 if the bathymetry is constant
-- the total modified energy
+- the total modified entropy/energy (integral of ``\hat U``/``\hat e``), which is called here
+  [`entropy_modified`](@ref) and [`energy_total_modified`](@ref)
 
 for periodic boundary conditions (see Ranocha and Ricchiuto (2025)).
 Additionally, it is well-balanced for the lake-at-rest stationary solution, see
@@ -784,7 +785,7 @@ 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 for flat bathymetry
-# - the total modified energy
+# - the total modified entropy/modified energy (integral of ``\hat U``/``\hat e``) as a nonlinear invariant
 # for periodic boundary conditions, see
 # - Hendrik Ranocha and Mario Ricchiuto (2025)
 #   Structure-Preserving Approximations of the Serre-Green-Naghdi
@@ -1352,22 +1353,22 @@ end
 
 # The modified entropy/energy takes the whole `q` for every point in space
 """
-    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SerreGreenNaghdiEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e_mod, q_global, equations::SerreGreenNaghdiEquations1D, cache)
 
-Return the modified total energy `e` of the primitive variables `q_global` for the
+Return the modified total energy ``\\hat e`` 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`
+non-hydrostatic contributions. The [`energy_total_modified`](@ref)
 is a conserved quantity (for periodic boundary conditions).
 
 For a [`bathymetry_flat`](@ref) the total modified energy is given by
 ```math
-\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2.
+\\hat e(\\eta, v) = \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2.
 ```
 For a [`bathymetry_mild_slope`](@ref) the total modified energy is given by
 ```math
-\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h (-h v_x + 1.5 v b_x)^2.
+\\hat e(\\eta, v) = \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h (-h v_x + 1.5 v b_x)^2.
 ```
 For a [`bathymetry_variable`](@ref) the total modified energy has the additional term
 ```math
@@ -1379,7 +1380,7 @@ For a [`bathymetry_variable`](@ref) the total modified energy has the additional
 
 See also [`energy_total_modified`](@ref).
 """
-function energy_total_modified!(e, q_global,
+function energy_total_modified!(e_mod, q_global,
                                 equations::SerreGreenNaghdiEquations1D,
                                 cache)
     # unpack physical parameters and SBP operator `D1`
@@ -1412,12 +1413,12 @@ function energy_total_modified!(e, q_global,
         # This is equivalent to the above expression when `integrate`d over
         # the domain, i.e., multiplied by 1^T M from the left.
         @.. D2.b = (1 / 6) * h^3
-        mul!(e, D2, v)
-        scale_by_mass_matrix!(e, D1)
-        @.. e = -e * v
-        scale_by_inverse_mass_matrix!(e, D1)
+        mul!(e_mod, D2, v)
+        scale_by_mass_matrix!(e_mod, D1)
+        @.. e_mod = -e_mod * v
+        scale_by_inverse_mass_matrix!(e_mod, D1)
 
-        @.. e = e + 1 / 2 * g * eta^2 + 1 / 2 * h * v^2
+        @.. e_mod = e_mod + 1 / 2 * g * eta^2 + 1 / 2 * h * v^2
     else
         # Periodic boundary conditions
         # and reflecting boundary conditions
@@ -1442,12 +1443,12 @@ function energy_total_modified!(e, q_global,
             end
         end
 
-        @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 +
-                1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2
+        @.. e_mod = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 +
+                    1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2
         if equations.bathymetry_type isa BathymetryVariable
-            @.. e += 1 / 8 * h * (v * b_x)^2
+            @.. e_mod += 1 / 8 * h * (v * b_x)^2
         end
     end
 
-    return e
+    return e_mod
 end
diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl
index 2d54c420..5f8d1264 100644
--- a/src/equations/svaerd_kalisch_1d.jl
+++ b/src/equations/svaerd_kalisch_1d.jl
@@ -33,7 +33,8 @@ The equations by Svärd and Kalisch are presented and analyzed in Svärd and Kal
 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 for flat bathymetry
-- the total modified energy
+- the total modified entropy/energy (integral of ``\hat U``/``\hat e``), which is called here
+  [`entropy_modified`](@ref) and [`energy_total_modified`](@ref)
 
 for periodic boundary conditions (see Lampert, Ranocha).
 Additionally, it is well-balanced for the lake-at-rest stationary solution, see Lampert and Ranocha (2025).
@@ -366,7 +367,7 @@ end
 # Discretization that conserves
 # - the total water (integral of ``h``) as a linear invariant
 # - the total momentum (integral of ``h v``) as a nonlinear invariant for flat bathymetry
-# - the total modified energy
+# - the total modified entropy/energy (integral of ``\hat U``/``\hat e``) as a nonlinear invariant
 # for periodic boundary conditions, see
 # - Joshua Lampert, Hendrik Ranocha (2025)
 #   Structure-preserving numerical methods for two nonlinear systems of dispersive wave equations
@@ -513,9 +514,9 @@ end
 
 # The modified entropy/energy takes the whole `q` for every point in space
 """
-    DispersiveShallowWater.energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D, cache)
+    DispersiveShallowWater.energy_total_modified!(e_mod, q_global, equations::SvaerdKalischEquations1D, cache)
 
-Return the modified total energy `e` of the primitive variables `q_global` for the
+Return the modified total energy ``\\hat e`` of the primitive variables `q_global` for the
 [`SvaerdKalischEquations1D`](@ref). It contains an additional term containing a
 derivative compared to the usual [`energy_total`](@ref) modeling
 non-hydrostatic contributions. The `energy_total_modified`
@@ -523,7 +524,7 @@ is a conserved quantity (for periodic boundary conditions).
 
 It is given by
 ```math
-\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta v_x^2.
+\\hat e(\\eta, v) = \\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{2} \\hat\\beta v_x^2.
 ```
 
 `q_global` is a vector of the primitive variables at ALL nodes.
@@ -531,8 +532,8 @@ It is given by
 
 See also [`energy_total_modified`](@ref).
 """
-@inline function energy_total_modified!(e, q_global, equations::SvaerdKalischEquations1D,
-                                        cache)
+@inline function energy_total_modified!(e_mod, q_global,
+                                        equations::SvaerdKalischEquations1D, cache)
     # unpack physical parameters and SBP operator `D1`
     g = gravity(equations)
     (; D1, h, b, v_x, beta_hat) = cache
@@ -550,6 +551,6 @@ See also [`energy_total_modified`](@ref).
         mul!(v_x, D1, v)
     end
 
-    @.. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 2 * beta_hat * v_x^2
-    return e
+    @.. e_mod = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 2 * beta_hat * v_x^2
+    return e_mod
 end