diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index f5d2f7b0bad..e900fd64073 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -301,6 +301,46 @@ end
     return abs(v1) + c, abs(v2) + c
 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::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_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 (we have p = kappa * rho^gamma)
+    c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1))
+    c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1))
+
+    λ_max = 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::PolytropicEulerEquations2D)
+    rho_ll, v1_ll, v2_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr = cons2prim(u_rr, equations)
+
+    # Calculate normal velocities and sound speed (we have p = kappa * rho^gamma)
+    # left
+    v_ll = (v1_ll * normal_direction[1] +
+            v2_ll * normal_direction[2])
+    c_ll = sqrt(equations.gamma * equations.kappa * rho_ll^(equations.gamma - 1))
+    # right
+    v_rr = (v1_rr * normal_direction[1] +
+            v2_rr * normal_direction[2])
+    c_rr = sqrt(equations.gamma * equations.kappa * rho_rr^(equations.gamma - 1))
+
+    return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
+end
+
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::PolytropicEulerEquations2D)
     rho, rho_v1, rho_v2 = u
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 3b8dc3c4331..c1379587cc8 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -858,6 +858,30 @@ end
     end
 end
 
+@timed_testset "Consistency check for Lax-Friedrich 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_lax_friedrichs(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_lax_friedrichs(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)