diff --git a/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl
new file mode 100644
index 00000000000..95ef38eb701
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_eulerpolytropic_convergence.jl
@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                periodicity = true,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.1)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index d4902bbafb2..f5d2f7b0bad 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -120,7 +120,7 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
     return prim2cons(SVector(rho, v1, v2), equations)
 end
 
-# Calculate 1D flux for a single point in the normal direction
+# Calculate 2D flux for a single point in the normal direction
 # Note, this directional vector is not normalized
 @inline function flux(u, normal_direction::AbstractVector,
                       equations::PolytropicEulerEquations2D)
@@ -135,8 +135,28 @@ end
     return SVector(f1, f2, f3)
 end
 
+# Calculate 2D flux for a single point
+@inline function flux(u, orientation::Integer, equations::PolytropicEulerEquations2D)
+    _, v1, v2 = cons2prim(u, equations)
+    p = pressure(u, equations)
+
+    rho_v1 = u[2]
+    rho_v2 = u[3]
+
+    if orientation == 1
+        f1 = rho_v1
+        f2 = rho_v1 * v1 + p
+        f3 = rho_v1 * v2
+    else
+        f1 = rho_v2
+        f2 = rho_v2 * v1
+        f3 = rho_v2 * v2 + p
+    end
+    return SVector(f1, f2, f3)
+end
+
 """
-    flux_winters_etal(u_ll, u_rr, normal_direction,
+    flux_winters_etal(u_ll, u_rr, orientation_or_normal_direction,
                       equations::PolytropicEulerEquations2D)
 
 Entropy conserving two-point flux for isothermal or polytropic gases.
@@ -178,6 +198,37 @@ For details see Section 3.2 of the following reference
     return SVector(f1, f2, f3)
 end
 
+@inline function flux_winters_etal(u_ll, u_rr, orientation::Integer,
+                                   equations::PolytropicEulerEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    # Compute the necessary mean values
+    if equations.gamma == 1.0 # isothermal gas
+        rho_mean = ln_mean(rho_ll, rho_rr)
+    else # equations.gamma > 1 # polytropic gas
+        rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
+    end
+    v1_avg = 0.5 * (v1_ll + v1_rr)
+    v2_avg = 0.5 * (v2_ll + v2_rr)
+    p_avg = 0.5 * (p_ll + p_rr)
+
+    if orientation == 1 # x-direction
+        f1 = rho_mean * 0.5 * (v1_ll + v1_rr)
+        f2 = f1 * v1_avg + p_avg
+        f3 = f1 * v2_avg
+    else # y-direction
+        f1 = rho_mean * 0.5 * (v2_ll + v2_rr)
+        f2 = f1 * v1_avg
+        f3 = f1 * v2_avg + p_avg
+    end
+
+    return SVector(f1, f2, f3)
+end
+
 @inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
                                      equations::PolytropicEulerEquations2D)
     rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
@@ -196,6 +247,53 @@ end
     return lambda_min, lambda_max
 end
 
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    if orientation == 1 # x-direction
+        λ_min = min(v1_ll - c_ll, v1_rr - c_rr)
+        λ_max = max(v1_ll + c_ll, v1_rr + c_rr)
+    else # y-direction
+        λ_min = min(v2_ll - c_ll, v2_rr - c_rr)
+        λ_max = max(v2_ll + c_ll, v2_rr + c_rr)
+    end
+
+    return λ_min, λ_max
+end
+
+# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes
+@inline function min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector,
+                                     equations::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+    # Pressure for polytropic Euler
+    p_ll = equations.kappa * rho_ll^equations.gamma
+    p_rr = equations.kappa * rho_rr^equations.gamma
+
+    norm_ = norm(normal_direction)
+
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_
+
+    v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # The v_normals are already scaled by the norm
+    λ_min = min(v_normal_ll - c_ll, v_normal_rr - c_rr)
+    λ_max = max(v_normal_ll + c_ll, v_normal_rr + c_rr)
+
+    return λ_min, λ_max
+end
+
 @inline function max_abs_speeds(u, equations::PolytropicEulerEquations2D)
     rho, v1, v2 = cons2prim(u, equations)
     c = sqrt(equations.gamma * equations.kappa * rho^(equations.gamma - 1))
diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index dd2248e10b2..17ba001d5ec 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -598,6 +598,29 @@ end
     end
 end
 
+@trixi_testset "elixir_eulerpolytropic_convergence.jl: HLL(Davis)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
+                        solver=DGSEM(polydeg = 3,
+                                     surface_flux = FluxHLL(min_max_speed_davis),
+                                     volume_integral = VolumeIntegralFluxDifferencing(volume_flux)),
+                        l2=[
+                            0.0016689832177644243, 0.0025920263793104445,
+                            0.003281074494629298,
+                        ],
+                        linf=[
+                            0.01099488320190023, 0.013309526619350365,
+                            0.02008032661117909,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_eulerpolytropic_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"),
                         l2=[
diff --git a/test/test_tree_2d_eulerpolytropic.jl b/test/test_tree_2d_eulerpolytropic.jl
new file mode 100644
index 00000000000..545cf7274ff
--- /dev/null
+++ b/test/test_tree_2d_eulerpolytropic.jl
@@ -0,0 +1,35 @@
+module TestExamples2DEulerMulticomponent
+
+using Test
+using Trixi
+
+include("test_trixi.jl")
+
+EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
+
+@testset "Polytropic Euler" begin
+#! format: noindent
+
+@trixi_testset "elixir_eulerpolytropic_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_eulerpolytropic_convergence.jl"),
+                        l2=[
+                            0.0016689832177626373, 0.0025920263793094526,
+                            0.003281074494626679,
+                        ],
+                        linf=[
+                            0.010994883201896677, 0.013309526619350365,
+                            0.02008032661117376,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+end
+
+end # module
diff --git a/test/test_tree_2d_part2.jl b/test/test_tree_2d_part2.jl
index d8e86d14f18..622f12109ff 100644
--- a/test/test_tree_2d_part2.jl
+++ b/test/test_tree_2d_part2.jl
@@ -26,6 +26,9 @@ isdir(outdir) && rm(outdir, recursive = true)
     # Compressible Euler Multicomponent
     include("test_tree_2d_eulermulti.jl")
 
+    # Compressible Polytropic Euler
+    include("test_tree_2d_eulerpolytropic.jl")
+
     # Compressible Euler coupled with acoustic perturbation equations
     include("test_tree_2d_euleracoustics.jl")
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index aa7840ec7aa..20368481607 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -782,6 +782,54 @@ end
     end
 end
 
+@timed_testset "Consistency check for HLL flux with Davis wave speed estimates: Polytropic CEE" begin
+    flux_hll = FluxHLL(min_max_speed_davis)
+
+    gamma = 1.4
+    kappa = 0.5     # Scaling factor for the pressure.
+    equations = PolytropicEulerEquations2D(gamma, kappa)
+    u = SVector(1.1, -0.5, 2.34)
+
+    orientations = [1, 2]
+    for orientation in orientations
+        @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+    end
+
+    normal_directions = [SVector(1.0, 0.0),
+        SVector(0.0, 1.0),
+        SVector(0.5, -0.5),
+        SVector(-1.2, 0.3)]
+
+    for normal_direction in normal_directions
+        @test flux_hll(u, u, normal_direction, equations) ≈
+              flux(u, normal_direction, equations)
+    end
+end
+
+@timed_testset "Consistency check for Winters flux: Polytropic CEE" begin
+    for gamma in [1.4, 1.0, 5 / 3]
+        kappa = 0.5     # Scaling factor for the pressure.
+        equations = PolytropicEulerEquations2D(gamma, kappa)
+        u = SVector(1.1, -0.5, 2.34)
+
+        orientations = [1, 2]
+        for orientation in orientations
+            @test flux_winters_etal(u, u, orientation, equations) ≈
+                  flux(u, orientation, equations)
+        end
+
+        normal_directions = [SVector(1.0, 0.0),
+            SVector(0.0, 1.0),
+            SVector(0.5, -0.5),
+            SVector(-1.2, 0.3)]
+
+        for normal_direction in normal_directions
+            @test flux_winters_etal(u, u, normal_direction, equations) ≈
+                  flux(u, normal_direction, equations)
+        end
+    end
+end
+
 @timed_testset "Consistency check for HLL flux with Davis wave speed estimates: LEE" begin
     flux_hll = FluxHLL(min_max_speed_davis)