diff --git a/examples/euler/dry_air/buoyancy/elixir_internal_energy_inertia_gravity_waves.jl b/examples/euler/dry_air/buoyancy/elixir_internal_energy_inertia_gravity_waves.jl
new file mode 100644
index 00000000..4709591b
--- /dev/null
+++ b/examples/euler/dry_air/buoyancy/elixir_internal_energy_inertia_gravity_waves.jl
@@ -0,0 +1,108 @@
+# 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::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    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 = CompressibleEulerInternalEnergyEquationsWithGravity2D(c_p = 1004,
+                                                                  c_v = 717,
+                                                                  gravity = 9.81)
+
+surface_flux = (flux_conservative_es, flux_nonconservative_es)
+volume_flux = (flux_conservative_etec, flux_nonconservative_etec)
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+boundary_conditions = (;
+                       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)
+trees_per_dimension = (60, 8)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = polydeg,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (true, false))
+
+initial_condition = initial_condition_gravity_waves
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    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 b81ae9f8..9d72ed0c 100644
--- a/src/TrixiAtmo.jl
+++ b/src/TrixiAtmo.jl
@@ -58,7 +58,8 @@ export CompressibleMoistEulerEquations2D,
        CompressibleEulerPotentialTemperatureEquationsWithGravity2D,
        CompressibleEulerPotentialTemperatureEquationsWithGravity3D,
        CompressibleEulerEnergyEquationsWithGravity2D,
-       CompressibleEulerEnergyEquationsWithGravity3D
+       CompressibleEulerEnergyEquationsWithGravity3D,
+       CompressibleEulerInternalEnergyEquationsWithGravity2D
 
 export GlobalCartesianCoordinates, GlobalSphericalCoordinates
 
@@ -72,7 +73,8 @@ export flux_nonconservative_zeros, flux_nonconservative_ec,
        flux_tec, flux_etec, flux_nonconservative_souza_etal,
        flux_nonconservative_artiano_etal,
        flux_nonconservative_waruszewski_etal, flux_zero,
-       flux_ec_rain, flux_LMARS
+       flux_ec_rain, flux_LMARS, flux_nonconservative_es, flux_conservative_es,
+       flux_conservative_etec, flux_nonconservative_etec
 
 export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier!
 
diff --git a/src/equations/compressible_euler_internal_energy_with_gravity_2d.jl b/src/equations/compressible_euler_internal_energy_with_gravity_2d.jl
new file mode 100644
index 00000000..1e927636
--- /dev/null
+++ b/src/equations/compressible_euler_internal_energy_with_gravity_2d.jl
@@ -0,0 +1,326 @@
+# 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""" 
+CompressibleEulerInternalEnergyEquationsWithGravity2D(; c_p, c_v, gravity) 
+
+The compressible Euler equations with gravity
+```math
+\frac{\partial}{\partial t}
+\begin{pmatrix}
+\rho \\ \rho v_1 \\ \rho v_2 \\ \rho e_{\text{internal}}
+\end{pmatrix}
++
+\frac{\partial}{\partial x}
+\begin{pmatrix}
+ \rho v_1 \\ \rho v_1^2 + p \\ \rho v_1 v_2 \\ (\rho e_{\text{internal}} +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_{\text{internal}} +p) v_2
+\end{pmatrix}
++
+\begin{pmatrix}
+0 \\ 0 \\ 0 \\ - v_1 \frac{\partial p}{\partial x} - v_2 \frac{\partial p}{\partial y}
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\ - \rho \frac{\partial}{\partial x} \phi \\ - \rho \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_{\text{internal}}`` is the specific internal energy, ``\phi`` is the gravitational potential, and
+```math
+p = (\gamma - 1) \rho e_{\text{internal}}
+```
+the pressure.
+"""
+struct CompressibleEulerInternalEnergyEquationsWithGravity2D{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 CompressibleEulerInternalEnergyEquationsWithGravity2D(; 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),
+                  ::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    ("rho", "rho_v1", "rho_v2", "rho_e_internal", "phi")
+end
+
+varnames(::typeof(cons2prim),
+::CompressibleEulerInternalEnergyEquationsWithGravity2D) = ("rho",
+                                                            "v1",
+                                                            "v2",
+                                                            "p", "phi")
+
+have_nonconservative_terms(::CompressibleEulerInternalEnergyEquationsWithGravity2D) = True()
+
+@inline function boundary_condition_slip_wall(u_inner,
+                                              normal_direction::AbstractVector,
+                                              x, t,
+                                              surface_flux_functions,
+                                              equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    # normalize the outward pointing direction
+    normal = normal_direction / norm(normal_direction)
+    surface_flux_function, nonconservative_flux_function = surface_flux_functions
+    # compute the normal momentum
+    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
+
+"""
+	flux_conservative_etec(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+
+Entropy conserving, total energy conserving and kinetic energy preserving two-point flux by
+-  Marco Artiano, Hendrik Ranocha (2026)
+   On Affordable High-Order Entropy-Conservative/Stable and 
+   Well-Balanced Methods for Nonconservative Hyperbolic Systems
+   [DOI: 10.48550/arXiv.2603.18978](https://arxiv.org/abs/2603.18978)
+"""
+@inline function flux_conservative_etec(u_ll, u_rr,
+                                        normal_direction::AbstractVector,
+                                        equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr, phi_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]
+    v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    # Compute the necessary mean values
+    rho_mean = ln_mean(rho_ll, rho_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
+
+    # 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 * inv_rho_p_mean * equations.inv_gamma_minus_one + p_avg * v_dot_n_avg
+    return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
+end
+
+"""
+	flux_nonconservative_etec(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+
+Nonconservative part of the entropy conserving, total energy conserving and kinetic energy preserving two-point flux by
+-  Marco Artiano, Hendrik Ranocha (2026)
+   On Affordable High-Order Entropy-Conservative/Stable and 
+   Well-Balanced Methods for Nonconservative Hyperbolic Systems
+  [DOI: 10.48550/arXiv.2603.18978](https://arxiv.org/abs/2603.18978)
+"""
+@inline function flux_nonconservative_etec(u_ll, u_rr,
+                                           normal_direction::AbstractVector,
+                                           equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr, phi_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]
+    v_dot_n_avg = 0.5f0 * (v_dot_n_rr + v_dot_n_ll)
+
+    # Compute the necessary mean values
+    rho_mean = ln_mean(rho_ll, rho_rr)
+
+    phi_jump = phi_rr - phi_ll
+    gravity = rho_mean * phi_jump
+    v_p_mean = -v_dot_n_avg * (p_rr - p_ll)
+    return SVector(zero(eltype(u_ll)), gravity * normal_direction[1],
+                   gravity * normal_direction[2], v_p_mean, zero(eltype(u_ll)))
+end
+
+"""
+	flux_conservative_es(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+
+Entropy stable two-point flux by
+-  Marco Artiano, Hendrik Ranocha (2026)
+   On Affordable High-Order Entropy-Conservative/Stable and 
+   Well-Balanced Methods for Nonconservative Hyperbolic Systems
+   [DOI: 10.48550/arXiv.2603.18978](https://arxiv.org/abs/2603.18978)
+"""
+@inline function flux_conservative_es(u_ll, u_rr, normal_direction::AbstractVector,
+                                      equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr, phi_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]
+    v_dot_n_avg = 0.5f0 * (v_dot_n_rr + v_dot_n_ll)
+
+    # Compute the necessary mean values
+    rho_mean = ln_mean(rho_ll, rho_rr)
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    a = 0.5f0 * (c_ll + c_rr)
+    norm_ = norm(normal_direction)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    rho_avg = 0.5f0 * (rho_ll + rho_rr)
+    p_interface = 0.5f0 * (p_ll + p_rr) -
+                  0.5f0 * a * rho_avg * (v_dot_n_rr - v_dot_n_ll) / norm_
+    v_interface = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) -
+                  1 / (2 * a * rho_avg) * (p_rr - p_ll) * norm_
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    if (v_interface >= 0)
+        rho_upwind = rho_mean - 0.5f0 * (rho_rr - rho_ll)
+        _, f2, f3, _, _ = u_ll * v_interface
+    else
+        rho_upwind = rho_mean + 0.5f0 * (rho_rr - rho_ll)
+        _, f2, f3, _, _ = u_rr * v_interface
+    end
+    inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_upwind * v_interface
+    f2 = f2 + p_interface * normal_direction[1]
+    f3 = f3 + p_interface * normal_direction[2]
+    f4 = f1 * inv_rho_p_mean * equations.inv_gamma_minus_one + p_avg * v_interface
+
+    return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
+end
+
+"""
+	flux_nonconservative_es(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+
+Nonconservative part of the entropy stable two-point flux by
+-  Marco Artiano, Hendrik Ranocha (2026)
+   On Affordable High-Order Entropy-Conservative/Stable and 
+   Well-Balanced Methods for Nonconservative Hyperbolic Systems
+   [DOI: 10.48550/arXiv.2603.18978](https://arxiv.org/abs/2603.18978)
+"""
+@inline function flux_nonconservative_es(u_ll, u_rr, normal_direction::AbstractVector,
+                                         equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll, phi_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr, phi_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]
+    v_dot_n_avg = 0.5f0 * (v_dot_n_rr + v_dot_n_ll)
+
+    # Compute the necessary mean values
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+    a = 0.5f0 * (c_ll + c_rr)
+    norm_ = norm(normal_direction)
+    rho_avg = 0.5f0 * (rho_ll + rho_rr)
+    v_interface = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) -
+                  1 / (2 * a * rho_avg) * (p_rr - p_ll) * norm_
+    p_avg = 0.5f0 * (p_ll + p_rr)
+    v_p_mean = -v_interface * (p_rr - p_ll)
+
+    return SVector(zero(eltype(u_ll)), zero(eltype(u_ll)), zero(eltype(u_ll)), v_p_mean,
+                   zero(eltype(u_ll)))
+end
+
+@inline function prim2cons(prim,
+                           equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    rho, v1, v2, p, phi = prim
+    rho_v1 = rho * v1
+    rho_v2 = rho * v2
+    rho_e_internal = p * equations.inv_gamma_minus_one
+    return SVector(rho, rho_v1, rho_v2, rho_e_internal, phi)
+end
+
+@inline function cons2prim(u,
+                           equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e_internal, phi = u
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * rho_e_internal
+    return SVector(rho, v1, v2, p, phi)
+end
+
+@inline function cons2cons(u,
+                           equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    return u
+end
+
+@inline function cons2entropy(u,
+                              equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    rho, rho_v1, rho_v2, rho_e_internal, phi = u
+
+    w1 = log(rho_e_internal * (equations.gamma - 1) / rho^equations.gamma) -
+         equations.gamma
+    w4 = rho / (rho_e_internal * (equations.gamma - 1))
+
+    return SVector(w1, zero(eltype(u)), zero(eltype(u)), w4, zero(eltype(u)))
+end
+
+@inline function entropy(cons,
+                         equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    p = pressure(cons, equations)
+    # Thermodynamic entropy
+    s = log(p) - equations.gamma * log(cons[1])
+    S = -s * cons[1] / (equations.gamma - 1)
+    return S
+end
+
+@inline function pressure(cons,
+                          equations::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    p = (equations.gamma - 1) * cons[4]
+    return p
+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::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    rho_ll, v1_ll, v2_ll, p_ll, phi = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr, phi = 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::CompressibleEulerInternalEnergyEquationsWithGravity2D)
+    rho, v1, v2, p = cons2prim(u, equations)
+    c = sqrt(equations.gamma * p / rho)
+
+    return abs(v1) + c, abs(v2) + c
+end
+end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index ee91577f..e045f68c 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -350,6 +350,7 @@ 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("compressible_euler_energy_with_gravity_3d.jl")
+include("compressible_euler_internal_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 08cb6874..4a3a57d9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -109,14 +109,20 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
     end
 
     @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_energy_2d"
-        @testset verbose=true showtiming=true "Euler internal energy 2D tests" begin
+        @testset verbose=true showtiming=true "Euler internal plus kinetic energy 2D tests" begin
             include("test_2d_euler_energy.jl")
         end
     end
 
     @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_energy_3d"
-        @testset verbose=true showtiming=true "Euler internal energy 3D tests" begin
+        @testset verbose=true showtiming=true "Euler internal plus kinetic energy 3D tests" begin
             include("test_3d_euler_energy.jl")
         end
     end
+
+    @time if TRIXI_TEST == "all" || TRIXI_TEST == "euler_internal energy_2d"
+        @testset verbose=true showtiming=true "Euler internal energy 2D tests" begin
+            include("test_2d_euler_internal_energy.jl")
+        end
+    end
 end
diff --git a/test/test_2d_euler_internal_energy.jl b/test/test_2d_euler_internal_energy.jl
new file mode 100644
index 00000000..c7c1b1b5
--- /dev/null
+++ b/test/test_2d_euler_internal_energy.jl
@@ -0,0 +1,29 @@
+module TestExamples2DEulerInternalEnergy
+
+include("test_trixiatmo.jl")
+
+EXAMPLES_DIR = joinpath(EXAMPLES_DIR, "euler/dry_air")
+
+@trixi_testset "elixir_internal_energy_inertia_gravity_waves" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "buoyancy",
+                                 "elixir_energy_inertia_gravity_waves.jl"),
+                        l2=[
+                            2.3801001547126282e-7,
+                            6.190703738154071e-6,
+                            3.382129115733546e-5,
+                            0.04382810096988029,
+                            9.251160452856857e-12
+                        ],
+                        linf=[
+                            1.2362858290426715e-6,
+                            6.0616987330064376e-5,
+                            0.00033470830317765867,
+                            0.27272386007825844,
+                            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 617ece0d..81bba1bd 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -580,4 +580,54 @@ end
                                              equations)) == RealT
     end
 end
+
+@timed_testset "Compressible Euler Internal Energy With Gravity 2D" begin
+    for RealT in (Float32, Float64)
+        equations = @inferred CompressibleEulerInternalEnergyEquationsWithGravity2D(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),
+                                                   RealT(2), one(RealT))
+        surface_flux_functions = (flux_conservative_etec,
+                                  flux_nonconservative_etec)
+        normal_direction = SVector(one(RealT), one(RealT))
+
+        @test eltype(@inferred flux_nonconservative_etec(u_ll, u_rr,
+                                                         normal_direction,
+                                                         equations)) ==
+              RealT
+        @test eltype(@inferred flux_conservative_etec(u_ll, u_rr,
+                                                      normal_direction,
+                                                      equations)) ==
+              RealT
+        @test eltype(@inferred flux_nonconservative_es(u_ll, u_rr,
+                                                       normal_direction,
+                                                       equations)) ==
+              RealT
+        @test eltype(@inferred flux_conservative_es(u_ll, u_rr,
+                                                    normal_direction,
+                                                    equations)) ==
+              RealT
+        @test eltype(@inferred boundary_condition_slip_wall(u_inner,
+                                                            normal_direction,
+                                                            x, t,
+                                                            surface_flux_functions,
+                                                            equations)) ==
+              SVector{5, RealT}
+        @test @inferred Trixi.have_nonconservative_terms(equations) === Trixi.True()
+        @test eltype(varnames(cons2cons, equations)) == String
+        @test eltype(varnames(cons2prim, equations)) == String
+        @test eltype(@inferred cons2prim(u, equations)) == RealT
+        @test eltype(@inferred cons2cons(u, equations)) == RealT
+        @test eltype(@inferred prim2cons(u, equations)) == RealT
+        @test eltype(@inferred cons2entropy(u, equations)) == RealT
+        @test typeof(@inferred pressure(u, equations)) == RealT
+        @test typeof(@inferred entropy(cons, equations)) == RealT
+        @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+        @test typeof(@inferred max_abs_speed(u_ll, u_rr, normal_direction,
+                                             equations)) == RealT
+    end
+end
 end