diff --git a/examples/elixir_euler_energy_inertia_gravity_waves.jl b/examples/elixir_euler_energy_inertia_gravity_waves.jl
new file mode 100644
index 00000000..304072a0
--- /dev/null
+++ b/examples/elixir_euler_energy_inertia_gravity_waves.jl
@@ -0,0 +1,112 @@
+# This test case is used to compute convergence rates via a linearized solution.
+# The setup follows the approach commonly adopted in benchmark studies; therefore,
+# a fixed CFL number is employed.
+#
+# References:
+# - Michael Baldauf and Slavko Brdar (2013):
+#   "An analytic solution for linear gravity waves in a channel as a test
+#   for numerical models using the non-hydrostatic, compressible Euler equations"
+#   Q. J. R. Meteorol. Soc., DOI: 10.1002/qj.2105
+#   https://doi.org/10.1002/qj.2105
+#
+# - Maciej Waruszewski, Jeremy E. Kozdon, Lucas C. Wilcox, Thomas H. Gibson,
+#   and Francis X. Giraldo (2022):
+#   "Entropy stable discontinuous Galerkin methods for balance laws
+#   in non-conservative form: Applications to the Euler equations with gravity"
+#   JCP, DOI: 10.1016/j.jcp.2022.111507
+#   https://doi.org/10.1016/j.jcp.2022.111507
+#
+# - Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025):
+#   "Structure-Preserving High-Order Methods for the Compressible Euler Equations
+#   in Potential Temperature Formulation for Atmospheric Flows"
+#   https://arxiv.org/abs/2509.10311
+
+using OrdinaryDiffEqSSPRK
+using Trixi, TrixiAtmo
+
+"""
+	initial_condition_gravity_waves(x, t,
+                                        equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Test cases for linearized analytical solution by
+-  Baldauf, Michael and Brdar, Slavko (2013)
+   An analytic solution for linear gravity waves in a channel as a test 
+   for numerical models using the non-hydrostatic, compressible {E}uler equations
+   [DOI: 10.1002/qj.2105] (https://doi.org/10.1002/qj.2105)
+"""
+function initial_condition_gravity_waves(x, t,
+                                         equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    g = equations.g
+    c_p = equations.c_p
+    c_v = equations.c_v
+    # center of perturbation
+    x_c = 100_000.0
+    a = 5_000
+    H = 10_000
+    R = c_p - c_v    # gas constant (dry air)
+    T0 = 250
+    delta = g / (R * T0)
+    DeltaT = 0.001
+    Tb = DeltaT * sinpi(x[2] / H) * exp(-(x[1] - x_c)^2 / a^2)
+    ps = 100_000  # reference pressure
+    rhos = ps / (T0 * R)
+    rho_b = rhos * (-Tb / T0)
+    p = ps * exp(-delta * x[2])
+    rho = rhos * exp(-delta * x[2]) + rho_b * exp(-0.5 * delta * x[2])
+    v1 = 20
+    v2 = 0
+    return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
+end
+
+equations = CompressibleEulerEnergyEquationsWithGravity2D(c_p = 1004,
+                                                          c_v = 717,
+                                                          gravity = 9.81)
+
+# We have an isothermal background state with T0 = 250 K. 
+# The reference speed of sound can be computed as:
+# cs = sqrt(gamma * R * T0)
+cs = sqrt(equations.gamma * equations.R * 250)
+surface_flux = (FluxLMARS(cs), flux_zero)
+volume_flux = (flux_ranocha, flux_nonconservative_waruzewski_etal)
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+boundary_conditions = (x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (300_000.0, 10_000.0)
+cells_per_dimension = (60, 8)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = (true, false))
+source_terms = nothing
+initial_condition = initial_condition_gravity_waves
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms,
+                                    boundary_conditions = boundary_conditions)
+tspan = (0.0, 1800.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 10000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        stepsize_callback)
+
+sol = solve(ode,
+            SSPRK43(thread = Trixi.True());
+            maxiters = 1.0e7,
+            dt = 1e-1, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks, adaptive = false)
diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl
index bdf1c846..1b7d8aa4 100644
--- a/src/TrixiAtmo.jl
+++ b/src/TrixiAtmo.jl
@@ -30,7 +30,8 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
                         energy_kinetic, energy_internal, energy_total, entropy, pressure,
                         flux, flux_ec, flux_chandrashekar, flux_wintermeyer_etal,
                         flux_fjordholm_etal, flux_nonconservative_wintermeyer_etal,
-                        flux_nonconservative_fjordholm_etal, FluxLMARS
+                        flux_nonconservative_fjordholm_etal, FluxLMARS, flux_shima_etal,
+                        flux_ranocha, flux_kennedy_gruber
 
 using Trixi: ln_mean, stolarsky_mean, inv_ln_mean
 
@@ -49,7 +50,8 @@ export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
        CompressibleEulerPotentialTemperatureEquations3D,
        CompressibleEulerPotentialTemperatureEquationsWithGravity1D,
        CompressibleEulerPotentialTemperatureEquationsWithGravity2D,
-       CompressibleEulerPotentialTemperatureEquationsWithGravity3D
+       CompressibleEulerPotentialTemperatureEquationsWithGravity3D,
+       CompressibleEulerEnergyEquationsWithGravity2D
 
 export GlobalCartesianCoordinates, GlobalSphericalCoordinates
 
diff --git a/src/equations/compressible_euler_energy_with_gravity_2d.jl b/src/equations/compressible_euler_energy_with_gravity_2d.jl
new file mode 100644
index 00000000..5419c83a
--- /dev/null
+++ b/src/equations/compressible_euler_energy_with_gravity_2d.jl
@@ -0,0 +1,644 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    CompressibleEulerEnergyEquationsWithGravity2D(gamma)
+
+The compressible Euler equations with gravity
+```math
+\frac{\partial}{\partial t}
+\begin{pmatrix}
+\rho \\ \rho v_1 \\ \rho v_2 \\ \rho e
+\end{pmatrix}
++
+\frac{\partial}{\partial x}
+\begin{pmatrix}
+ \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e +p) v_1
+\end{pmatrix}
++
+\frac{\partial}{\partial y}
+\begin{pmatrix}
+\rho v_2 \\ \rho v_1 v_2 \\ \rho v_2^2 + p \\ (\rho e +p) v_2
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\ - rho \frac{\partial}{\partial x} \phi \\ - rho \frac{\partial}{\partial y} \phi \\ - rho v \frac{\partial}{\partial x} \phi - rho v \frac{\partial}{\partial y} \phi
+\end{pmatrix}
+```
+for an ideal gas with ratio of specific heats `gamma`
+in two space dimensions.
+Here, ``\rho`` is the density, ``v_1``, ``v_2`` are the velocities, ``e`` is the specific total energy **rather than** specific internal energy, ``phi`` is the gravitational potential, and
+```math
+p = (\gamma - 1) \left( \rho e - \frac{1}{2} \rho (v_1^2+v_2^2) \right)
+```
+the pressure.
+"""
+struct CompressibleEulerEnergyEquationsWithGravity2D{RealT <: Real} <:
+       AbstractCompressibleEulerEquations{2, 5}
+    p_0::RealT # reference pressure in Pa
+    c_p::RealT # specific heat at constant pressure in J/(kg K)
+    c_v::RealT # specific heat at constant volume in J/(kg K)
+    g::RealT # gravitational acceleration in m/s²
+    R::RealT # gas constant in J/(kg K)
+    gamma::RealT # ratio of specific heats 
+    inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
+    function CompressibleEulerEnergyEquationsWithGravity2D(; c_p, c_v,
+                                                           gravity)
+        c_p, c_v, g = promote(c_p, c_v, gravity)
+        p_0 = 100_000
+        R = c_p - c_v
+        gamma = c_p / c_v
+        inv_gamma_minus_one = inv(gamma - 1)
+        return new{typeof(c_p)}(p_0, c_p, c_v, g, R,
+                                gamma,
+                                inv_gamma_minus_one)
+    end
+end
+
+function varnames(::typeof(cons2cons),
+                  ::CompressibleEulerEnergyEquationsWithGravity2D)
+    ("rho", "rho_v1", "rho_v2", "rho_e", "phi")
+end
+varnames(::typeof(cons2prim), ::CompressibleEulerEnergyEquationsWithGravity2D) = ("rho",
+                                                                                  "v1",
+                                                                                  "v2",
+                                                                                  "p",
+                                                                                  "phi")
+
+have_nonconservative_terms(::CompressibleEulerEnergyEquationsWithGravity2D) = Trixi.True()
+
+# Slip-wall boundary condition
+# Determine the boundary numerical surface flux for a slip wall condition.
+# Imposes a zero normal velocity at the wall.
+@inline function boundary_condition_slip_wall(u_inner,
+                                              normal_direction::AbstractVector,
+                                              x, t,
+                                              surface_flux_functions,
+                                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # normalize the outward pointing direction
+    normal = normal_direction / norm(normal_direction)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    # compute the normal velocity
+    u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3]
+
+    # create the "external" boundary solution state
+    u_boundary = SVector(u_inner[1],
+                         u_inner[2] - 2 * u_normal * normal[1],
+                         u_inner[3] - 2 * u_normal * normal[2],
+                         u_inner[4], u_inner[5])
+
+    # calculate the boundary flux
+    flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+    noncons_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
+                                                 equations)
+    return flux, noncons_flux
+end
+
+@inline function boundary_condition_slip_wall(u_inner, orientation,
+                                              direction, x, t,
+                                              surface_flux_functions,
+                                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+
+    ## get the appropriate normal vector from the orientation
+    if orientation == 1
+        u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4],
+                             u_inner[5])
+    else # orientation == 2
+        u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4],
+                             u_inner[5])
+    end
+
+    # Calculate boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_inner, u_boundary, orientation,
+                                                     equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+        noncons_flux = nonconservative_flux_function(u_boundary, u_inner, orientation,
+                                                     equations)
+    end
+
+    return flux, noncons_flux
+end
+
+# Calculate 1D flux for a single point
+@inline function flux(u, orientation::Integer,
+                      equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e = u
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2))
+    if orientation == 1
+        f1 = rho_v1
+        f2 = rho_v1 * v1 + p
+        f3 = rho_v1 * v2
+        f4 = (rho_e + p) * v1
+    else
+        f1 = rho_v2
+        f2 = rho_v2 * v1
+        f3 = rho_v2 * v2 + p
+        f4 = (rho_e + p) * v2
+    end
+    return SVector(f1, f2, f3, f4, zero(eltype(u)))
+end
+
+# Calculate 1D flux for a single point in the normal direction
+# Note, this directional vector is not normalized
+@inline function flux(u, normal_direction::AbstractVector,
+                      equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_e = last(u)
+    rho, v1, v2, p = cons2prim(u, equations)
+
+    v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
+    rho_v_normal = rho * v_normal
+    f1 = rho_v_normal
+    f2 = rho_v_normal * v1 + p * normal_direction[1]
+    f3 = rho_v_normal * v2 + p * normal_direction[2]
+    f4 = (rho_e + p) * v_normal
+    return SVector(f1, f2, f3, f4, zero(eltype(u)))
+end
+
+"""
+	flux_nonconservative_waruzewski_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Well-balanced gravity term for an isothermal background state
+for the [`CompressibleEulerEnergyEquationsWithGravity2D`](@ref)
+developed by
+
+-  Maciej Waruszewski and Jeremy E. Kozdon and Lucas C. Wilcox and Thomas H. Gibson and Francis X. Giraldo (2022)
+   Entropy stable discontinuous {G}alerkin methods for balance laws 
+   in non-conservative form: Applications to the {E}uler equations with gravity
+   [DOI: 10.1016/j.jcp.2022.111507](https://doi.org/10.1016/j.jcp.2022.111507)
+
+The well-balancedness on curvilinear coordinates was proven by
+-  Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025)
+   Structure-Preserving High-Order Methods for the Compressible Euler Equations 
+   in Potential Temperature Formulation for Atmospheric Flows
+   (https://arxiv.org/abs/2509.10311)
+"""
+@inline function flux_nonconservative_waruzewski_etal(u_ll, u_rr,
+                                                      normal_direction::AbstractVector,
+                                                      equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _, phi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _, phi_rr = u_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_ll = rho_v2_ll / rho_ll
+    v2_rr = rho_v2_rr / rho_rr
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    rho_avg = ln_mean(rho_ll, rho_rr)
+    phi_jump = phi_rr - phi_ll
+    return SVector(zero(eltype(u_ll)),
+                   normal_direction[1] * rho_avg * phi_jump,
+                   normal_direction[2] * rho_avg * phi_jump,
+                   normal_direction[1] * rho_avg * v1_avg * phi_jump +
+                   normal_direction[2] * rho_avg * v2_avg * phi_jump,
+                   zero(eltype(u_ll)))
+end
+
+"""
+	flux_nonconservative_artiano_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Well-balanced gravity term for a constant potential temperature background state
+for the [`CompressibleEulerEnergyEquationsWithGravity2D`](@ref)
+developed by
+
+-  Marco Artiano, Oswald Knoth, Peter Spichtinger, Hendrik Ranocha (2025)
+   Structure-Preserving High-Order Methods for the Compressible Euler Equations 
+   in Potential Temperature Formulation for Atmospheric Flows
+   (https://arxiv.org/abs/2509.10311)
+"""
+@inline function flux_nonconservative_artiano_etal(u_ll, u_rr,
+                                                   normal_direction::AbstractVector,
+                                                   equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _, phi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _, phi_rr = u_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_ll = rho_v2_ll / rho_ll
+    v2_rr = rho_v2_rr / rho_rr
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    rho_avg = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
+    phi_jump = phi_rr - phi_ll
+    return SVector(zero(eltype(u_ll)),
+                   normal_direction[1] * rho_avg * phi_jump,
+                   normal_direction[2] * rho_avg * phi_jump,
+                   normal_direction[1] * rho_avg * v1_avg * phi_jump +
+                   normal_direction[2] * rho_avg * v2_avg * phi_jump,
+                   zero(eltype(u_ll)))
+end
+
+"""
+	flux_nonconservative_souza_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Kinetic and potential energy preserving (KPEP) gravity term 
+for the [`CompressibleEulerEnergyEquationsWithGravity2D`](@ref)
+developed by
+
+-  Souza et al. 
+   The Flux-Differencing Discontinuous {G}alerkin Method Applied to 
+   an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere
+   [DOI: 10.1029/2022MS003527] (https://doi.org/10.1029/2022MS003527)
+"""
+@inline function flux_nonconservative_souza_etal(u_ll, u_rr,
+                                                 normal_direction::AbstractVector,
+                                                 equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, _, phi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, _, phi_rr = u_rr
+    v1_ll = rho_v1_ll / rho_ll
+    v1_rr = rho_v1_rr / rho_rr
+    v2_ll = rho_v2_ll / rho_ll
+    v2_rr = rho_v2_rr / rho_rr
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    rho_avg = 0.5f0 * (rho_ll + rho_rr)
+    phi_jump = phi_rr - phi_ll
+    return SVector(zero(eltype(u_ll)),
+                   normal_direction[1] * rho_avg * phi_jump,
+                   normal_direction[2] * rho_avg * phi_jump,
+                   normal_direction[1] * rho_avg * v1_avg * phi_jump +
+                   normal_direction[2] * rho_avg * v2_avg * phi_jump,
+                   zero(eltype(u_ll)))
+end
+
+"""
+    flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction,
+                    equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+This flux is is a modification of the original kinetic energy preserving two-point flux by
+- Yuichi Kuya, Kosuke Totani and Soshi Kawai (2018)
+  Kinetic energy and entropy preserving schemes for compressible flows
+  by split convective forms
+  [DOI: 10.1016/j.jcp.2018.08.058](https://doi.org/10.1016/j.jcp.2018.08.058)
+
+The modification is in the energy flux to guarantee pressure equilibrium and was developed by
+- Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, Soshi Kawai (JCP 2020)
+  Preventing spurious pressure oscillations in split convective form discretizations for
+  compressible flows
+  [DOI: 10.1016/j.jcp.2020.110060](https://doi.org/10.1016/j.jcp.2020.110060)
+"""
+@inline function flux_shima_etal(u_ll, u_rr, normal_direction::AbstractVector,
+                                 equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # Average each factor of products in flux
+    rho_avg = 0.5f0 * (rho_ll + rho_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_avg * v_dot_n_avg
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+    f4 = (f1 * velocity_square_avg +
+          p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
+          + 0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
+
+    return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
+end
+
+"""
+    flux_kennedy_gruber(u_ll, u_rr, orientation_or_normal_direction,
+                        equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Kinetic energy preserving two-point flux by
+- Kennedy and Gruber (2008)
+  Reduced aliasing formulations of the convective terms within the
+  Navier-Stokes equations for a compressible fluid
+  [DOI: 10.1016/j.jcp.2007.09.020](https://doi.org/10.1016/j.jcp.2007.09.020)
+"""
+@inline function flux_kennedy_gruber(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_e_ll = last(u_ll)
+    rho_e_rr = last(u_rr)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Average each factor of products in flux
+    rho_avg = 0.5f0 * (rho_ll + rho_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    v_dot_n_avg = v1_avg * normal_direction[1] + v2_avg * normal_direction[2]
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    e_avg = 0.5f0 * (rho_e_ll / rho_ll + rho_e_rr / rho_rr)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_avg * v_dot_n_avg
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+    f4 = f1 * e_avg + p_avg * v_dot_n_avg
+
+    return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
+end
+
+"""
+    flux_ranocha(u_ll, u_rr, orientation_or_normal_direction,
+                 equations::CompressibleEulerEnergyEquationsWithGravity2D)
+
+Entropy conserving and kinetic energy preserving two-point flux by
+- Hendrik Ranocha (2018)
+  Generalised Summation-by-Parts Operators and Entropy Stability of Numerical Methods
+  for Hyperbolic Balance Laws
+  [PhD thesis, TU Braunschweig](https://cuvillier.de/en/shop/publications/7743)
+See also
+- Hendrik Ranocha (2020)
+  Entropy Conserving and Kinetic Energy Preserving Numerical Methods for
+  the Euler Equations Using Summation-by-Parts Operators
+  [Proceedings of ICOSAHOM 2018](https://doi.org/10.1007/978-3-030-39647-3_42)
+"""
+@inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector,
+                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # Compute the necessary mean values
+    rho_mean = ln_mean(rho_ll, rho_rr)
+    # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
+    # in exact arithmetic since
+    #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
+    #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
+    inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+    f4 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one)
+          +
+          0.5f0 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
+
+    return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
+end
+
+@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector,
+                                         equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    c = flux_lmars.speed_of_sound
+
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # Note that this is the same as computing v_ll and v_rr with a normalized normal vector
+    # and then multiplying v by `norm_` again, but this version is slightly faster.
+    norm_ = norm(normal_direction)
+
+    rho = 0.5f0 * (rho_ll + rho_rr)
+    p = 0.5f0 * (p_ll + p_rr) - 0.5f0 * c * rho * (v_rr - v_ll) / norm_
+    v = 0.5f0 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
+
+    # We treat the energy term analogous to the potential temperature term in the paper by
+    # Chen et al., i.e. we use p_ll and p_rr, and not p
+    if v >= 0
+        f1, f2, f3, f4 = u_ll * v
+        f4 = f4 + p_ll * v
+    else
+        f1, f2, f3, f4 = u_rr * v
+        f4 = f4 + p_rr * v
+    end
+
+    return SVector(f1,
+                   f2 + p * normal_direction[1],
+                   f3 + p * normal_direction[2],
+                   f4, zero(eltype(u_ll)))
+end
+
+# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
+# maximum velocity magnitude plus the maximum speed of sound
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Get the velocity value in the appropriate direction
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    else # orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    end
+    # Calculate sound speeds
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
+end
+
+@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speed
+    # left
+    v_ll = (v1_ll * normal_direction[1]
+            +
+            v2_ll * normal_direction[2])
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    # right
+    v_rr = (v1_rr * normal_direction[1]
+            +
+            v2_rr * normal_direction[2])
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
+end
+
+# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
+                               equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Get the velocity value in the appropriate direction
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    else # orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    end
+    # Calculate sound speeds
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    return max(abs(v_ll) + c_ll, abs(v_rr) + c_rr)
+end
+
+# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
+@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
+                               equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speeds
+    # left
+    v_ll = (v1_ll * normal_direction[1]
+            +
+            v2_ll * normal_direction[2])
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    # right
+    v_rr = (v1_rr * normal_direction[1]
+            +
+            v2_rr * normal_direction[2])
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    norm_ = norm(normal_direction)
+    return max(abs(v_ll) + c_ll * norm_,
+               abs(v_rr) + c_rr * norm_)
+end
+
+@inline function max_abs_speeds(u,
+                                equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, v1, v2, p = cons2prim(u, equations)
+    c = sqrt(equations.gamma * p / rho)
+
+    return abs(v1) + c, abs(v2) + c
+end
+
+# Convert conservative variables to primitive
+@inline function cons2prim(u,
+                           equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e, phi = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1 * v1 + rho_v2 * v2))
+
+    return SVector(rho, v1, v2, p, phi)
+end
+
+# Convert conservative variables to entropy
+@inline function cons2entropy(u,
+                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v_square = v1^2 + v2^2
+    p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho * v_square)
+    s = log(p) - equations.gamma * log(rho)
+    rho_p = rho / p
+
+    w1 = (equations.gamma - s) * equations.inv_gamma_minus_one -
+         0.5f0 * rho_p * v_square
+    w2 = rho_p * v1
+    w3 = rho_p * v2
+    w4 = -rho_p
+
+    return SVector(w1, w2, w3, w4, zero(eltype(u)))
+end
+
+@inline function entropy2cons(w,
+                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD
+    # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1)
+    @unpack gamma = equations
+
+    # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986)
+    # instead of `-rho * s / (gamma - 1)`
+    V1, V2, V3, V5 = w .* (gamma - 1)
+
+    # s = specific entropy, eq. (53)
+    s = gamma - V1 + (V2^2 + V3^2) / (2 * V5)
+
+    # eq. (52)
+    rho_iota = ((gamma - 1) / (-V5)^gamma)^(equations.inv_gamma_minus_one) *
+               exp(-s * equations.inv_gamma_minus_one)
+
+    # eq. (51)
+    rho = -rho_iota * V5
+    rho_v1 = rho_iota * V2
+    rho_v2 = rho_iota * V3
+    rho_e = rho_iota * (1 - (V2^2 + V3^2) / (2 * V5))
+    return SVector(rho, rho_v1, rho_v2, rho_e, zero(eltype(w)))
+end
+
+# Convert primitive to conservative variables
+@inline function prim2cons(prim,
+                           equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, v1, v2, p, phi = prim
+    rho_v1 = rho * v1
+    rho_v2 = rho * v2
+    rho_e = p * equations.inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2)
+    return SVector(rho, rho_v1, rho_v2, rho_e, phi)
+end
+
+@inline function pressure(u,
+                          equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e = u
+    p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho)
+    return p
+end
+
+# Calculate thermodynamic entropy for a conservative state `cons`
+@inline function entropy_thermodynamic(cons,
+                                       equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # Pressure
+    p = pressure(cons, equations)
+
+    # Thermodynamic entropy
+    s = log(p) - equations.gamma * log(cons[1])
+
+    return s
+end
+
+# Calculate mathematical entropy for a conservative state `cons`
+@inline function entropy_math(cons,
+                              equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    # Mathematical entropy
+    S = -entropy_thermodynamic(cons, equations) * cons[1] *
+        equations.inv_gamma_minus_one
+
+    return S
+end
+
+# Default entropy is the mathematical entropy
+@inline function entropy(cons,
+                         equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    entropy_math(cons, equations)
+end
+
+# Calculate total energy for a conservative state `cons`
+@inline energy_total(cons, ::CompressibleEulerEnergyEquationsWithGravity2D) = cons[4]
+
+# Calculate kinetic energy for a conservative state `cons`
+@inline function energy_kinetic(u,
+                                equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e = u
+    return (rho_v1^2 + rho_v2^2) / (2 * rho)
+end
+
+# Calculate internal energy for a conservative state `cons`
+@inline function energy_internal(cons,
+                                 equations::CompressibleEulerEnergyEquationsWithGravity2D)
+    return energy_total(cons, equations) - energy_kinetic(cons, equations)
+end
+end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 151341fb..3e48ae3a 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -344,6 +344,7 @@ include("compressible_euler_potential_temperature_3d.jl")
 include("compressible_euler_potential_temperature_gravity_1d.jl")
 include("compressible_euler_potential_temperature_gravity_2d.jl")
 include("compressible_euler_potential_temperature_gravity_3d.jl")
+include("compressible_euler_energy_with_gravity_2d.jl")
 include("shallow_water_3d.jl")
 include("reference_data.jl")
 end # @muladd
diff --git a/test/runtests.jl b/test/runtests.jl
index c1447abc..a481633a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -70,6 +70,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
         include("test_3d_potential_temperature.jl")
     end
 
+    @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_energy_2d"
+        include("test_2d_euler_energy.jl")
+    end
+
     @time if TRIXI_TEST == "upstream"
         include("test_trixi_consistency.jl")
     end
diff --git a/test/test_2d_euler_energy.jl b/test/test_2d_euler_energy.jl
new file mode 100644
index 00000000..069f5599
--- /dev/null
+++ b/test/test_2d_euler_energy.jl
@@ -0,0 +1,27 @@
+module TestExamples2DEulerEnergy
+
+include("test_trixiatmo.jl")
+
+@trixi_testset "elixir_euler_energy_inertia_gravity_waves_2d" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_energy_inertia_gravity_waves.jl"),
+                        l2=[
+                            2.3800999105272615e-7,
+                            6.190703408721927e-6,
+                            3.3821288686132935e-5,
+                            0.04382810097184301,
+                            9.251160452856857e-12
+                        ],
+                        linf=[
+                            1.2362858294867607e-6,
+                            6.0616987290984525e-5,
+                            0.00033470830315131473,
+                            0.27272386016556993,
+                            4.3655745685100555e-11
+                        ], tspan=(0.0, 10.0), atol=5e-11)
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100)
+end
+
+end
diff --git a/test/test_type.jl b/test/test_type.jl
index 90402c69..e70881cd 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -368,5 +368,79 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                                        equations)) == RealT
         end
     end
+
+    @timed_testset "Compressible Euler Energy With Gravity 2D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred CompressibleEulerEnergyEquationsWithGravity2D(c_p = RealT(1004),
+                                                                                c_v = RealT(717),
+                                                                                gravity = RealT(9.81))
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT),
+                                                       one(RealT), one(RealT))
+
+            normal_direction = SVector(one(RealT), one(RealT))
+            surface_flux_function = (flux_lax_friedrichs, flux_zero)
+            orientations = [1, 2]
+            directions = [1, 2]
+
+            for direction in directions, orientation in orientations
+                @test eltype(@inferred boundary_condition_slip_wall(u_inner, orientation,
+                                                                    direction,
+                                                                    x, t,
+                                                                    surface_flux_function,
+                                                                    equations)) ==
+                      SVector{5, RealT}
+                @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation,
+                                                           equations)) == RealT
+                @test typeof(@inferred max_abs_speed(u_ll, u_rr, orientation,
+                                                     equations)) == RealT
+            end
+            @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
+            @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, normal_direction,
+                                                       equations)) ==
+                  RealT
+            @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_artiano_etal(u_ll, u_rr,
+                                                                     normal_direction,
+                                                                     equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_waruzewski_etal(u_ll, u_rr,
+                                                                        normal_direction,
+                                                                        equations)) ==
+                  RealT
+            @test eltype(@inferred flux_nonconservative_souza_etal(u_ll, u_rr,
+                                                                   normal_direction,
+                                                                   equations)) ==
+                  RealT
+            @test eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                                normal_direction,
+                                                                x, t,
+                                                                surface_flux_function,
+                                                                equations)) ==
+                  SVector{5, RealT}
+            @test eltype(varnames(cons2prim, equations)) == String
+
+            @test eltype(@inferred cons2prim(u, equations)) == RealT
+            @test eltype(@inferred prim2cons(u, equations)) == RealT
+            @test eltype(@inferred cons2entropy(u, equations)) == RealT
+            @test eltype(@inferred entropy2cons(cons2entropy(u, equations), equations)) ==
+                  RealT
+            @test typeof(@inferred pressure(u, equations)) == RealT
+            @test typeof(@inferred entropy(cons, equations)) == RealT
+            @test typeof(@inferred energy_kinetic(cons, equations)) == RealT
+            @test typeof(@inferred energy_total(cons, equations)) == RealT
+            @test typeof(@inferred energy_internal(cons, equations)) == RealT
+            @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+            @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
+                                                       equations)) == RealT
+            @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction,
+                                                 equations)) == RealT
+        end
+    end
 end
 end