diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
new file mode 100644
index 00000000000..cc4535be028
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -0,0 +1,85 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+"""
+    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+
+A discontinuous initial condition taken from
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)   
+"""
+function initial_condition_discontinuity(x, t,
+                                         equations::CompressibleEulerEquationsQuasi1D)
+    rho = (x[1] < 0) ? 3.4718 : 2.0
+    v1 = (x[1] < 0) ? -2.5923 : -3.0
+    p = (x[1] < 0) ? 5.7118 : 2.639
+    a = (x[1] < 0) ? 1.0 : 1.5
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_discontinuity
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+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.5)
+
+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/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
new file mode 100644
index 00000000000..ae1b2b24b62
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -0,0 +1,73 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations with a discontinuous nozzle width function.
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+# Setup a truly discontinuous density function and nozzle width for
+# this academic testcase of entropy conservation. The errors from the analysis
+# callback are not important but the entropy error for this test case
+# `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=6`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D)
+    v1 = 0.1
+    rho = 2.0 + 0.1 * x[1]
+    p = 3.0
+    a = 2.0 + x[1]
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_ec
+
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+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.8)
+
+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/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
new file mode 100644
index 00000000000..91bb1ba6e8c
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -0,0 +1,60 @@
+using OrdinaryDiffEq
+using Trixi
+using ForwardDiff
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = -1.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+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.8)
+
+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/Trixi.jl b/src/Trixi.jl
index b8110cf5bdd..e7b849e2642 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -139,6 +139,7 @@ export AcousticPerturbationEquations2D,
        CompressibleEulerEquations3D,
        CompressibleEulerMulticomponentEquations1D,
        CompressibleEulerMulticomponentEquations2D,
+       CompressibleEulerEquationsQuasi1D,
        IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
        IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D,
diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
new file mode 100644
index 00000000000..0a543277ee4
--- /dev/null
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -0,0 +1,328 @@
+# 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"""
+    CompressibleEulerEquationsQuasi1D(gamma)
+
+The quasi-1d compressible Euler equations (see Chan et al.  [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)  for details)
+```math
+\frac{\partial}{\partial t}
+\begin{pmatrix}
+a \rho \\ a \rho v_1 \\ a e
+\end{pmatrix}
++
+\frac{\partial}{\partial x}
+\begin{pmatrix}
+a \rho v_1 \\ a \rho v_1^2 \\ a v_1 (e +p)
+\end{pmatrix}
++ 
+a \frac{\partial}{\partial x}
+\begin{pmatrix}
+0 \\ p \\ 0    
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\ 0 \\ 0
+\end{pmatrix}
+```
+for an ideal gas with ratio of specific heats `gamma` in one space dimension.
+Here, ``\rho`` is the density, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, 
+``a`` the (possibly) variable nozzle width, and
+```math
+p = (\gamma - 1) \left( e - \frac{1}{2} \rho v_1^2 \right)
+```
+the pressure.
+
+The nozzle width function ``a(x)`` is set inside the initial condition routine
+for a particular problem setup. To test the conservative form of the compressible Euler equations one can set the 
+nozzle width variable ``a`` to one. 
+
+In addition to the unknowns, Trixi.jl currently stores the nozzle width values at the approximation points 
+despite being fixed in time.
+This affects the implementation and use of these equations in various ways:
+* The flux values corresponding to the nozzle width must be zero.
+* The nozzle width values must be included when defining initial conditions, boundary conditions or
+  source terms.
+* [`AnalysisCallback`](@ref) analyzes this variable.
+* Trixi.jl's visualization tools will visualize the nozzle width by default.
+"""
+struct CompressibleEulerEquationsQuasi1D{RealT <: Real} <:
+       AbstractCompressibleEulerEquations{1, 4}
+    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 CompressibleEulerEquationsQuasi1D(gamma)
+        γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1))
+        new{typeof(γ)}(γ, inv_gamma_minus_one)
+    end
+end
+
+have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = True()
+function varnames(::typeof(cons2cons), ::CompressibleEulerEquationsQuasi1D)
+    ("a_rho", "a_rho_v1", "a_e", "a")
+end
+function varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D)
+    ("rho", "v1", "p", "a")
+end
+
+"""
+    initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
+
+A smooth initial condition used for convergence tests in combination with
+[`source_terms_convergence_test`](@ref)
+(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
+"""
+function initial_condition_convergence_test(x, t,
+                                            equations::CompressibleEulerEquationsQuasi1D)
+    c = 2
+    A = 0.1
+    L = 2
+    f = 1 / L
+    ω = 2 * pi * f
+    ini = c + A * sin(ω * (x[1] - t))
+
+    rho = ini
+    v1 = 1.0
+    e = ini^2 / rho
+    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
+    a = 1.5 - 0.5 * cos(x[1] * pi)
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+"""
+    source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquationsQuasi1D)
+
+Source terms used for convergence tests in combination with
+[`initial_condition_convergence_test`](@ref)
+(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
+
+This manufactured solution source term is specifically designed for the mozzle width 'a(x) = 1.5 - 0.5 * cos(x[1] * pi)'
+as defined in [`initial_condition_convergence_test`](@ref).
+"""
+@inline function source_terms_convergence_test(u, x, t,
+                                               equations::CompressibleEulerEquationsQuasi1D)
+    # Same settings as in `initial_condition_convergence_test`. 
+    # Derivatives calculated with ForwardDiff.jl
+    c = 2
+    A = 0.1
+    L = 2
+    f = 1 / L
+    ω = 2 * pi * f
+    x1, = x
+    ini(x1, t) = c + A * sin(ω * (x1 - t))
+
+    rho(x1, t) = ini(x1, t)
+    v1(x1, t) = 1.0
+    e(x1, t) = ini(x1, t)^2 / rho(x1, t)
+    p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2)
+    a(x1, t) = 1.5 - 0.5 * cos(x1 * pi)
+
+    arho(x1, t) = a(x1, t) * rho(x1, t)
+    arhou(x1, t) = arho(x1, t) * v1(x1, t)
+    aE(x1, t) = a(x1, t) * e(x1, t)
+
+    darho_dt(x1, t) = ForwardDiff.derivative(t -> arho(x1, t), t)
+    darhou_dx(x1, t) = ForwardDiff.derivative(x1 -> arhou(x1, t), x1)
+
+    arhouu(x1, t) = arhou(x1, t) * v1(x1, t)
+    darhou_dt(x1, t) = ForwardDiff.derivative(t -> arhou(x1, t), t)
+    darhouu_dx(x1, t) = ForwardDiff.derivative(x1 -> arhouu(x1, t), x1)
+    dp1_dx(x1, t) = ForwardDiff.derivative(x1 -> p1(x1, t), x1)
+
+    auEp(x1, t) = a(x1, t) * v1(x1, t) * (e(x1, t) + p1(x1, t))
+    daE_dt(x1, t) = ForwardDiff.derivative(t -> aE(x1, t), t)
+    dauEp_dx(x1, t) = ForwardDiff.derivative(x1 -> auEp(x1, t), x1)
+
+    du1 = darho_dt(x1, t) + darhou_dx(x1, t)
+    du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t)
+    du3 = daE_dt(x1, t) + dauEp_dx(x1, t)
+
+    return SVector(du1, du2, du3, 0.0)
+end
+
+# Calculate 1D flux for a single point
+@inline function flux(u, orientation::Integer,
+                      equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho, v1, p, a = cons2prim(u, equations)
+    e = a_e / a
+
+    # Ignore orientation since it is always "1" in 1D
+    f1 = a_rho_v1
+    f2 = a_rho_v1 * v1
+    f3 = a * v1 * (e + p)
+
+    return SVector(f1, f2, f3, zero(eltype(u)))
+end
+
+"""
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+
+Non-symmetric two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the pressure  [`CompressibleEulerEquationsQuasi1D`](@ref) 
+and the nozzle width.
+
+Further details are available in the paper:
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)    
+"""
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+    #Variables
+    _, _, p_ll, a_ll = cons2prim(u_ll, equations)
+    _, _, p_rr, _ = cons2prim(u_rr, equations)
+
+    # For flux differencing using non-conservative terms, we return the 
+    # non-conservative flux scaled by 2. This cancels with a factor of 0.5 
+    # in the arithmetic average of {p}.
+    p_avg = p_ll + p_rr
+
+    z = zero(eltype(u_ll))
+
+    return SVector(z, a_ll * p_avg, z, z)
+end
+
+"""
+@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
+                                           equations::CompressibleEulerEquationsQuasi1D)
+
+Conservative (symmetric) part of the entropy conservative flux for quasi 1D compressible Euler equations split form.
+This flux is a generalization of [`flux_ranocha`](@ref) for [`CompressibleEulerEquations1D`](@ref).
+Further details are available in the paper:
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) 
+  High order entropy stable schemes for the quasi-one-dimensional
+  shallow water and compressible Euler equations
+  [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)     
+"""
+@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
+                                equations::CompressibleEulerEquationsQuasi1D)
+    # Unpack left and right state
+    rho_ll, v1_ll, p_ll, a_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, p_rr, a_rr = cons2prim(u_rr, equations)
+
+    # 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.5 * (v1_ll + v1_rr)
+    a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr)
+    p_avg = 0.5 * (p_ll + p_rr)
+    velocity_square_avg = 0.5 * (v1_ll * v1_rr)
+
+    # Calculate fluxes
+    # Ignore orientation since it is always "1" in 1D
+    f1 = rho_mean * a_v1_avg
+    f2 = rho_mean * a_v1_avg * v1_avg
+    f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) +
+         0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll)
+
+    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+end
+
+# Calculate estimates for 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::CompressibleEulerEquationsQuasi1D)
+    a_rho_ll, a_rho_v1_ll, a_e_ll, a_ll = u_ll
+    a_rho_rr, a_rho_v1_rr, a_e_rr, a_rr = u_rr
+
+    # Calculate primitive variables and speed of sound
+    rho_ll = a_rho_ll / a_ll
+    e_ll = a_e_ll / a_ll
+    v1_ll = a_rho_v1_ll / a_rho_ll
+    v_mag_ll = abs(v1_ll)
+    p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2)
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    rho_rr = a_rho_rr / a_rr
+    e_rr = a_e_rr / a_rr
+    v1_rr = a_rho_v1_rr / a_rho_rr
+    v_mag_rr = abs(v1_rr)
+    p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
+end
+
+@inline function max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho = a_rho / a
+    v1 = a_rho_v1 / a_rho
+    e = a_e / a
+    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
+    c = sqrt(equations.gamma * p / rho)
+
+    return (abs(v1) + c,)
+end
+
+# Convert conservative variables to primitive. We use the convention that the primitive
+# variables for the quasi-1D equations are `(rho, v1, p)` (i.e., the same as the primitive
+# variables for `CompressibleEulerEquations1D`)
+@inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    q = cons2prim(SVector(a_rho, a_rho_v1, a_e) / a,
+                  CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(q[1], q[2], q[3], a)
+end
+
+# The entropy for the quasi-1D compressible Euler equations is the entropy for the
+# 1D compressible Euler equations scaled by the channel width `a`.
+@inline function entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    q = a * entropy(SVector(a_rho, a_rho_v1, a_e) / a,
+                CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(q[1], q[2], q[3], a)
+end
+
+# Convert conservative variables to entropy. The entropy variables for the 
+# quasi-1D compressible Euler equations are identical to the entropy variables
+# for the standard Euler equations for an appropriate definition of `entropy`.
+@inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    w = cons2entropy(SVector(a_rho, a_rho_v1, a_e) / a,
+                     CompressibleEulerEquations1D(equations.gamma))
+
+    # we follow the convention for other spatially-varying equations such as
+    # `ShallowWaterEquations1D` and return the spatially varying coefficient 
+    # `a` as the final entropy variable.
+    return SVector(w[1], w[2], w[3], a)
+end
+
+# Convert primitive to conservative variables
+@inline function prim2cons(u, equations::CompressibleEulerEquationsQuasi1D)
+    rho, v1, p, a = u
+    q = prim2cons(u, CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(a * q[1], a * q[2], a * q[3], a)
+end
+
+@inline function density(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, _, _, a = u
+    rho = a_rho / a
+    return rho
+end
+
+@inline function pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    return pressure(SVector(a_rho, a_rho_v1, a_e) / a,
+                    CompressibleEulerEquations1D(equations.gamma))
+end
+
+@inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a,
+                            CompressibleEulerEquations1D(equations.gamma))
+end
+end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 582d672b756..7a3c326984d 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -409,6 +409,7 @@ abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <:
 include("compressible_euler_1d.jl")
 include("compressible_euler_2d.jl")
 include("compressible_euler_3d.jl")
+include("compressible_euler_quasi_1d.jl")
 
 # CompressibleEulerMulticomponentEquations
 abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <:
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 6cd7998ab02..39a1f6e30ba 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -393,6 +393,79 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_quasi_1d_source_terms.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_source_terms.jl"),
+                        l2=[
+                            3.876288369618363e-7,
+                            2.2247043122302947e-7,
+                            2.964004224572679e-7,
+                            5.2716983399807875e-8,
+                        ],
+                        linf=[
+                            2.3925118561862746e-6,
+                            1.3603693522767912e-6,
+                            1.821888865105592e-6,
+                            1.1166012159335992e-7,
+                        ])
+    # 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_euler_quasi_1d_discontinuous.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_quasi_1d_discontinuous.jl"),
+                        l2=[
+                            0.045510421156346015,
+                            0.036750584788912195,
+                            0.2468985959132176,
+                            0.03684494180829024,
+                        ],
+                        linf=[
+                            0.3313374853025697,
+                            0.11621933362158643,
+                            1.827403013568638,
+                            0.28045939999015723,
+                        ])
+    # 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_euler_quasi_1d_ec.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_ec.jl"),
+                        l2=[
+                            0.08889113985713998,
+                            0.16199235348889673,
+                            0.40316524365054346,
+                            2.9602775074723667e-16,
+                        ],
+                        linf=[
+                            0.28891355898284043,
+                            0.3752709888964313,
+                            0.84477102402413,
+                            8.881784197001252e-16,
+                        ])
+    # 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_unit.jl b/test/test_unit.jl
index d2e744da62f..b3ed29d38e3 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -654,6 +654,19 @@ end
     end
 end
 
+@timed_testset "Consistency check for flux_chan_etal: CEEQ" begin
+
+    # Set up equations and dummy conservative variables state
+    equations = CompressibleEulerEquationsQuasi1D(1.4)
+    u = SVector(1.1, 2.34, 5.5, 2.73)
+
+    orientations = [1]
+    for orientation in orientations
+        @test flux_chan_etal(u, u, orientation, equations) ≈
+              flux(u, orientation, equations)
+    end
+end
+
 @timed_testset "Consistency check for HLL flux (naive): LEE" begin
     flux_hll = FluxHLL(min_max_speed_naive)