diff --git a/examples/tree_2d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl b/examples/tree_2d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl
new file mode 100644
index 00000000000..4dc4a3aa06a
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_density_wave_adaptive_vol_int.jl
@@ -0,0 +1,77 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations2D(1.4)
+
+# We repeat test case for linear stability of EC fluxes from the paper
+# - Gregor J. Gassner, Magnus Svärd, Florian J. Hindenlang (2020)
+#   Stability issues of entropy-stable and/or split-form high-order schemes
+#   [DOI: 10.1007/s10915-021-01720-8](https://doi.org/10.1007/s10915-021-01720-8)
+initial_condition = initial_condition_density_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_chandrashekar
+
+polydeg = 5
+basis = LobattoLegendreBasis(polydeg)
+
+volume_integral_weakform = VolumeIntegralWeakForm()
+volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(volume_flux)
+
+# This indicator compares the entropy production of the weak form to the 
+# true entropy evolution in that cell.
+# If the weak form dissipates more entropy than the true evolution
+# the indicator renders this admissible. Otherwise, the more stable
+# volume integral is to be used.
+indicator = IndicatorEntropyChange(maximum_entropy_increase = 0.0)
+
+# Adaptive volume integral using the entropy change indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(indicator = indicator,
+                                         volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff)
+
+#volume_integral = volume_integral_weakform # Stable, but unphysical entropy increase!
+#volume_integral = volume_integral_fluxdiff # Crashes!
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2, # 4 x 4 elements
+                n_cells_max = 30_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+# In the paper, CFL = 0.05 is used. Same observations recovered for CFL = 0.9, though.
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl
new file mode 100644
index 00000000000..081db374b8c
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl
@@ -0,0 +1,96 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    RealT = eltype(x)
+    slope = 15
+    B = tanh(slope * x[2] + 7.5f0) - tanh(slope * x[2] - 7.5f0)
+    rho = 0.5f0 + 0.75f0 * B
+    v1 = 0.5f0 * (B - 1)
+    v2 = convert(RealT, 0.1) * sinpi(2 * x[1])
+    p = 1
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+volume_integral_weakform = VolumeIntegralWeakForm()
+volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(volume_flux)
+
+# This indicator compares the entropy production of the weak form to the 
+# true entropy evolution in that cell.
+# If the weak form dissipates more entropy than the true evolution
+# the indicator renders this admissible. Otherwise, the more stable
+# volume integral is to be used.
+indicator = IndicatorEntropyChange(maximum_entropy_increase = 0.0)
+
+# Adaptive volume integral using the entropy change indicator to perform the 
+# stabilized/EC volume integral when needed and keeping the weak form if it is more diffusive.
+volume_integral = VolumeIntegralAdaptive(indicator = indicator,
+                                         volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff)
+
+#volume_integral = volume_integral_weakform # Crashes
+#volume_integral = volume_integral_fluxdiff # Crashes as well!
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 100_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.25)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_errors = Symbol[],
+                                     extra_analysis_integrals = (entropy,),
+                                     save_analysis = true)
+
+alive_callback = AliveCallback(alive_interval = 200)
+
+stepsize_callback = StepsizeCallback(cfl = 1.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 5e-3, ode_default_options()..., callback = callbacks);
diff --git a/src/Trixi.jl b/src/Trixi.jl
index baa391a0e4a..29105bd5406 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -251,6 +251,7 @@ export density, pressure, density_pressure, velocity, temperature,
        equilibrium_distribution,
        waterheight, waterheight_pressure
 export entropy, entropy_thermodynamic, entropy_math, entropy_guermond_etal,
+       entropy_potential,
        energy_total, energy_kinetic, energy_internal, energy_internal_specific,
        energy_magnetic, cross_helicity, magnetic_field, divergence_cleaning_field,
        enstrophy, vorticity
@@ -268,6 +269,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, VolumeIntegralShockCapturingRRG,
+       VolumeIntegralAdaptive, IndicatorEntropyChange,
        IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
index 0f7b2933c91..0f19ac2b02f 100644
--- a/src/callbacks_step/analysis_dg2d.jl
+++ b/src/callbacks_step/analysis_dg2d.jl
@@ -185,6 +185,68 @@ function calc_error_norms(func, u, t, analyzer,
     return l2_error, linf_error
 end
 
+# Use quadrature to numerically integrate a single element.
+# We do not multiply by the Jacobian to stay in reference space.
+# This avoids the need to divide the RHS of the DG scheme by the Jacobian when computing
+# the time derivative of entropy, see `entropy_change_reference_element`.
+function integrate_reference_element(func::Func, u, element,
+                                     mesh::AbstractMesh{2}, equations, dg::DGSEM, cache,
+                                     args...) where {Func}
+    @unpack weights = dg.basis
+
+    # Initialize integral with zeros of the right shape
+    element_integral = zero(func(u, 1, 1, element, equations, dg, args...))
+
+    for j in eachnode(dg), i in eachnode(dg)
+        element_integral += weights[i] * weights[j] *
+                            func(u, i, j, element, equations, dg, args...)
+    end
+
+    return element_integral
+end
+
+# Calculate ∫_e (∂S/∂u ⋅ ∂u/∂t) dΩ_e where the result on element 'e' is kept in reference space
+# Note that ∂S/∂u = w(u) with entropy variables w
+function entropy_change_reference_element(du::AbstractArray{<:Any, 4}, u,
+                                          element,
+                                          mesh::AbstractMesh{2},
+                                          equations, dg::DGSEM, cache, args...)
+    return integrate_reference_element(u, element, mesh, equations, dg, cache,
+                                       du) do u, i, j, element, equations, dg, du
+        u_node = get_node_vars(u, equations, dg, i, j, element)
+        du_node = get_node_vars(du, equations, dg, i, j, element)
+
+        dot(cons2entropy(u_node, equations), du_node)
+    end
+end
+
+# calculate surface integral of func(u, equations) * normal on the reference element.
+function surface_integral(func::Func, u, element,
+                          mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
+                          args...) where {Func}
+    @unpack weights = dg.basis
+
+    u_tmp = get_node_vars(u, equations, dg, 1, 1, element)
+    surface_integral = zero(func(u_tmp, 1, equations))
+    for i in eachnode(dg)
+        # integrate along x direction, normal in y (2) direction
+        u_bottom = get_node_vars(u, equations, dg, i, 1, element)
+        u_top = get_node_vars(u, equations, dg, i, nnodes(dg), element)
+
+        surface_integral += weights[i] *
+                            (func(u_top, 2, equations) - func(u_bottom, 2, equations))
+
+        # integrate along y direction, normal in x (1) direction
+        u_left = get_node_vars(u, equations, dg, 1, i, element)
+        u_right = get_node_vars(u, equations, dg, nnodes(dg), i, element)
+
+        surface_integral += weights[i] *
+                            (func(u_right, 1, equations) - func(u_left, 1, equations))
+    end
+
+    return surface_integral
+end
+
 function integrate_via_indices(func::Func, u,
                                mesh::TreeMesh{2}, equations, dg::DGSEM, cache,
                                args...; normalize = true) where {Func}
@@ -214,8 +276,8 @@ function integrate_via_indices(func::Func, u,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
                                            UnstructuredMesh2D, P4estMesh{2},
                                            T8codeMesh{2}},
-                               equations,
-                               dg::DGSEM, cache, args...; normalize = true) where {Func}
+                               equations, dg::DGSEM, cache,
+                               args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
     # Initialize integral with zeros of the right shape
diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 088e04588f3..cf733e97c3a 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -2235,6 +2235,28 @@ end
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
 
+@doc raw"""
+    entropy_potential(u, orientation::Int, 
+                      equations::AbstractCompressibleEulerEquations)
+
+Calculate the entropy potential, which for the compressible Euler equations is simply 
+the momentum for the choice of mathematical [`entropy`](@ref) ``S(u) = \frac{\rho s}{\gamma - 1}``
+with thermodynamic entropy ``s = \ln(p) - \gamma \ln(\rho)``.
+    
+## References
+- Eitan Tadmor (2003)
+  Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems
+  [DOI: 10.1017/S0962492902000156](https://doi.org/10.1017/S0962492902000156)
+"""
+@inline function entropy_potential(u, orientation::Int,
+                                   equations::CompressibleEulerEquations2D)
+    if orientation == 1
+        return u[2]
+    else # if orientation == 2
+        return u[3]
+    end
+end
+
 # State validation for Newton-bisection method of subcell IDP limiting
 @inline function Base.isvalid(u, equations::CompressibleEulerEquations2D)
     if u[1] <= 0 || pressure(u, equations) <= 0
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index 73d909567f8..ac7910374eb 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -282,6 +282,80 @@ function Base.show(io::IO, mime::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralAdaptive(;
+                           indicator = IndicatorEntropyChange(),
+                           volume_integral_default,
+                           volume_integral_stabilized)
+
+!!! warning "Experimental code"
+    This code is experimental and may change in any future release.
+
+The `indicator` is currently limited to [`IndicatorEntropyChange`](@ref).
+"""
+struct VolumeIntegralAdaptive{Indicator,
+                              VolumeIntegralDefault, VolumeIntegralStabilized} <:
+       AbstractVolumeIntegral
+    indicator::Indicator # A-posteriori indicator called after computation of `volume_integral_default`
+    volume_integral_default::VolumeIntegralDefault # Cheap(er) default volume integral to be used in non-critical regions
+    volume_integral_stabilized::VolumeIntegralStabilized # More expensive volume integral with stabilizing effect
+end
+
+function VolumeIntegralAdaptive(;
+                                indicator = IndicatorEntropyChange(),
+                                volume_integral_default,
+                                volume_integral_stabilized)
+    if !(indicator isa IndicatorEntropyChange)
+        throw(ArgumentError("`indicator` must be of type `IndicatorEntropyChange`."))
+    end
+
+    return VolumeIntegralAdaptive{typeof(indicator),
+                                  typeof(volume_integral_default),
+                                  typeof(volume_integral_stabilized)}(indicator,
+                                                                      volume_integral_default,
+                                                                      volume_integral_stabilized)
+end
+
+function Base.show(io::IO, mime::MIME"text/plain",
+                   integral::VolumeIntegralAdaptive)
+    @nospecialize integral # reduce precompilation time
+    @unpack volume_integral_default, volume_integral_stabilized, indicator = integral
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        summary_header(io, "VolumeIntegralAdaptive")
+
+        summary_line(io, "volume integral default",
+                     volume_integral_default |> typeof |> nameof)
+        if !(volume_integral_default isa VolumeIntegralWeakForm)
+            show(increment_indent(io), mime, volume_integral_default)
+        end
+
+        summary_line(io, "volume integral stabilized",
+                     volume_integral_stabilized |> typeof |> nameof)
+        show(increment_indent(io), mime, volume_integral_stabilized)
+
+        summary_line(io, "indicator", indicator |> typeof |> nameof)
+        show(increment_indent(io), mime, indicator)
+
+        summary_footer(io)
+    end
+end
+
+function create_cache(mesh, equations,
+                      volume_integral::VolumeIntegralAdaptive,
+                      dg, cache_containers, uEltype)
+    cache_default = create_cache(mesh, equations,
+                                 volume_integral.volume_integral_default,
+                                 dg, cache_containers, uEltype)
+    cache_stabilized = create_cache(mesh, equations,
+                                    volume_integral.volume_integral_stabilized,
+                                    dg, cache_containers, uEltype)
+
+    return (; cache_default..., cache_stabilized...)
+end
+
 # Abstract supertype for first-order `VolumeIntegralPureLGLFiniteVolume` and
 # second-order `VolumeIntegralPureLGLFiniteVolumeO2` subcell-based finite volume
 # volume integrals.
diff --git a/src/solvers/dgsem/calc_volume_integral.jl b/src/solvers/dgsem/calc_volume_integral.jl
index 7d71c97d051..e72ca1bf624 100644
--- a/src/solvers/dgsem/calc_volume_integral.jl
+++ b/src/solvers/dgsem/calc_volume_integral.jl
@@ -60,6 +60,43 @@ end
     return nothing
 end
 
+@inline function volume_integral_kernel!(du, u, element, mesh,
+                                         have_nonconservative_terms, equations,
+                                         volume_integral::VolumeIntegralAdaptive{<:IndicatorEntropyChange},
+                                         dg::DGSEM, cache)
+    @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral
+    @unpack maximum_entropy_increase = indicator
+
+    volume_integral_kernel!(du, u, element, mesh,
+                            have_nonconservative_terms, equations,
+                            volume_integral_default, dg, cache)
+
+    # Compute entropy production of the default volume integral.
+    # Minus sign because of the flipped sign of the volume term in the DG RHS.
+    # No scaling by inverse Jacobian here, as there is no Jacobian multiplication
+    # in `integrate_reference_element`.
+    dS_default = -entropy_change_reference_element(du, u, element,
+                                                   mesh, equations, dg, cache)
+
+    # Compute true entropy change given by surface integral of the entropy potential
+    dS_true = surface_integral(entropy_potential, u, element,
+                               mesh, equations, dg, cache)
+
+    entropy_change = dS_default - dS_true
+    if entropy_change > maximum_entropy_increase # Recompute using EC FD volume integral
+        # Reset default volume integral contribution.
+        # Note that this assumes that the volume terms are computed first,
+        # before any surface terms are added.
+        du[.., element] .= zero(eltype(du))
+
+        volume_integral_kernel!(du, u, element, mesh,
+                                have_nonconservative_terms, equations,
+                                volume_integral_stabilized, dg, cache)
+    end
+
+    return nothing
+end
+
 function calc_volume_integral!(du, u, mesh,
                                have_nonconservative_terms, equations,
                                volume_integral, dg::DGSEM, cache)
diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl
index 52d1a671030..3b706e0f09c 100644
--- a/src/solvers/dgsem/dgsem.jl
+++ b/src/solvers/dgsem/dgsem.jl
@@ -133,5 +133,6 @@ end
 end
 
 include("containers.jl")
+include("indicators.jl")
 include("calc_volume_integral.jl")
 end # @muladd
diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem/indicators.jl
similarity index 73%
rename from src/solvers/dgsem_tree/indicators.jl
rename to src/solvers/dgsem/indicators.jl
index fdc2ca3a536..9d7e7d6c514 100644
--- a/src/solvers/dgsem_tree/indicators.jl
+++ b/src/solvers/dgsem/indicators.jl
@@ -5,6 +5,8 @@
 @muladd begin
 #! format: noindent
 
+# Abstract supertype of indicators used for AMR, shock capturing, and
+# adaptive volume-integral selection
 abstract type AbstractIndicator end
 
 function create_cache(typ::Type{IndicatorType},
@@ -272,4 +274,93 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax)
         summary_box(io, "IndicatorMax", setup)
     end
 end
+
+@doc raw"""
+    IndicatorEntropyChange(; maximum_entropy_increase::Real = 0.0)
+
+This indicator checks the difference in mathematical [`entropy`](@ref) (``S``) due to the application
+of a volume integral (VI) compared to the true/analytical entropy evolution
+(without any dissipation inside the element).
+In particular, the indicator computes
+```math
+\Delta S = \dot{S}_\mathrm{VI} - \dot{S}_\text{true} =
+\int_{\Omega_m}
+\frac{\partial S}{\partial \boldsymbol{u}} \cdot \dot{\boldsymbol u}_\mathrm{VI}
+\mathrm{d} \Omega_m
+- 
+\int_{\partial \Omega_m}
+\boldsymbol{\psi} \cdot \hat{\boldsymbol{n}}
+\mathrm{d} \partial \Omega_m
+```
+for the currently processed element/cell ``m``.
+Here, ``\dot{\boldsymbol u}_\mathrm{VI}`` is the change in the DG right-hand-side due to the volume integral only.
+``\dot{S}_\text{true}`` is the true entropy evolution, which can be computed from the
+entropy potential ``\boldsymbol{\psi}`` (see also [`entropy_potential`](@ref)).
+
+This is discussed in more detail in
+- Chen, Shu (2017)
+  "Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws"
+  [DOI: 10.1016/j.jcp.2017.05.025](https://doi.org/10.1016/j.jcp.2017.05.025)
+- Lin, Chan (2024)
+  "High order entropy stable discontinuous Galerkin spectral element methods through subcell limiting"
+  [DOI: 10.1016/j.jcp.2023.112677](https://doi.org/10.1016/j.jcp.2023.112677)
+
+For ``\Delta S < \sigma \leq 0`` with ``\sigma`` being set to `maximum_entropy_increase`,
+the e.g. [`VolumeIntegralWeakForm`](@ref) is more entropy-diffusive than the true entropy change
+(which could be recovered with the [`VolumeIntegralFluxDifferencing`](@ref) and an
+entropy-conservative flux such as [`flux_ranocha`](@ref), for instance).
+
+If ``\sigma > 0`` is set, i.e., `maximum_entropy_increase > 0`, the indicator allows for
+limited entropy increase, thereby allowing to use e.g. the cheaper weak-form volume integral
+even in slightly entropy-producing situations to reduce computational cost.
+
+Supposed to be used in conjunction with [`VolumeIntegralAdaptive`](@ref) which then selects
+a volume integral for every cell/element ``m``.
+
+The logic behind this indicator is similar to the "companion" scheme
+approach proposed in Chapter 5 of
+
+- Carpenter, Fisher, Nielsen, and Frankel (2014)
+  "Entropy Stable Spectral Collocation Schemes for the Navier-Stokes Equations: Discontinuous Interfaces"
+  [DOI: 10.1137/130932193](https://doi.org/10.1137/130932193)
+
+Here, we thus equip e.g. the flux-differencing volume integral with a "companion" weak-form
+volume integral.
+However, usage of the entropy potential allows for comparison with the true entropy change.
+
+!!! note
+    This indicator is **not implemented as an AMR indicator**, i.e., it is **not
+    possible** to employ this as the `indicator` in [`ControllerThreeLevel`](@ref),
+    for instance.
+"""
+struct IndicatorEntropyChange{RealT <: Real} <:
+       AbstractIndicator
+    maximum_entropy_increase::RealT
+
+    function IndicatorEntropyChange(; maximum_entropy_increase = 0.0)
+        return new{typeof(maximum_entropy_increase)}(maximum_entropy_increase)
+    end
+end
+
+function Base.show(io::IO, indicator::IndicatorEntropyChange)
+    @nospecialize indicator # reduce precompilation time
+
+    print(io, "IndicatorEntropyChange(")
+    print(io, "maximum_entropy_increase=", indicator.maximum_entropy_increase, ")")
+
+    return nothing
+end
+
+function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorEntropyChange)
+    @nospecialize indicator # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, indicator)
+    else
+        setup = [
+            "maximum_entropy_increase" => indicator.maximum_entropy_increase
+        ]
+        summary_box(io, "IndicatorEntropyChange", setup)
+    end
+end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl
index af4615726b0..54d71af2b66 100644
--- a/src/solvers/dgsem_tree/dg.jl
+++ b/src/solvers/dgsem_tree/dg.jl
@@ -24,7 +24,6 @@ function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeri
 end
 
 # Indicators used for shock-capturing and AMR
-include("indicators.jl")
 include("indicators_1d.jl")
 include("indicators_2d.jl")
 include("indicators_3d.jl")
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 37df79487f5..f2b8a234278 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -575,6 +575,38 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_euler_density_wave_adaptive_vol_int.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_density_wave_adaptive_vol_int.jl"),
+                        l2=[
+                            0.034677655030611376,
+                            0.00346776550306035,
+                            0.0069355310061214445,
+                            0.0008669413757631413
+                        ],
+                        linf=[
+                            0.13358170277079373,
+                            0.013358170277035189,
+                            0.026716340554125667,
+                            0.0033395425693640846
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+
+    # Test/cover `show`
+    @test_nowarn show(stdout, indicator)
+    @test_nowarn show(IOContext(IOBuffer(), :compact => true), MIME"text/plain"(),
+                      indicator)
+    @test_nowarn show(IOContext(IOBuffer(), :compact => false), MIME"text/plain"(),
+                      indicator)
+
+    @test_nowarn show(IOContext(IOBuffer(), :compact => true), MIME"text/plain"(),
+                      volume_integral)
+    @test_nowarn show(IOContext(IOBuffer(), :compact => false), MIME"text/plain"(),
+                      volume_integral)
+end
+
 @trixi_testset "elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_kelvin_helmholtz_instability_fjordholm_etal.jl"),
@@ -669,6 +701,27 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 15000)
 end
 
+@trixi_testset "elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl"),
+                        tspan=(0.0, 0.1),
+                        l2=[
+                            0.026078588204610408,
+                            0.020356972101548874,
+                            0.02851024015440993,
+                            0.02951575049839344
+                        ],
+                        linf=[
+                            0.120658689916814,
+                            0.10672357779516586,
+                            0.06256826531635475,
+                            0.1199900961506577
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_euler_colliding_flow.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"),
                         l2=[