diff --git a/NEWS.md b/NEWS.md
index 3a3a504a911..02a723fca45 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -11,6 +11,7 @@ for human readability.
 - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`
 - Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, 
   can now be digested by Trixi in 2D and 3D.
+- Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh`
 
 ## Changes when updating to v0.6 from v0.5.x
 
diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl
new file mode 100644
index 00000000000..1817672778a
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl
@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+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
+    slope = 15
+    amplitude = 0.02
+    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+    rho = 0.5 + 0.75 * B
+    v1 = 0.5 * (B - 1)
+    v2 = 0.1 * sin(2 * pi * x[1])
+    p = 1.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+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 = 5,
+                n_cells_max = 100_000)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.7)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+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)
+
+save_restart = SaveRestartCallback(interval = 1000,
+                                   save_final_restart = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback,
+                        save_restart, save_solution)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl
index fe9ad92467f..74d0370647a 100644
--- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl
+++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl
@@ -22,7 +22,7 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
     r = sqrt(x[1]^2 + x[2]^2)
 
     pmax = 10.0
-    pmin = 1.0
+    pmin = 0.01
     rhomax = 1.0
     rhomin = 0.01
     if r <= 0.09
@@ -52,7 +52,8 @@ basis = LobattoLegendreBasis(3)
 
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 positivity_variables_cons = ["rho"],
-                                positivity_correction_factor = 0.5)
+                                positivity_variables_nonlinear = [pressure],
+                                positivity_correction_factor = 0.1)
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
@@ -84,7 +85,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-cfl = 0.5
+cfl = 0.4
 stepsize_callback = StepsizeCallback(cfl = cfl)
 
 glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl
index 9f34a6b3b4b..4dbf44d29c4 100644
--- a/src/callbacks_stage/subcell_bounds_check.jl
+++ b/src/callbacks_stage/subcell_bounds_check.jl
@@ -97,6 +97,9 @@ function init_callback(callback::BoundsCheckCallback, semi, limiter::SubcellLimi
                 end
                 print(f, ", " * string(variables[v]) * "_min")
             end
+            for variable in limiter.positivity_variables_nonlinear
+                print(f, ", " * string(variable) * "_min")
+            end
         end
         println(f)
     end
@@ -142,6 +145,11 @@ end
             println(string(variables[v]) * ":\n- positivity: ",
                     idp_bounds_delta_global[Symbol(string(v), "_min")])
         end
+        for variable in limiter.positivity_variables_nonlinear
+            variable_string = string(variable)
+            println(variable_string * ":\n- positivity: ",
+                    idp_bounds_delta_global[Symbol(variable_string, "_min")])
+        end
     end
     println("─"^100 * "\n")
 
diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl
index 545d19b5136..19d73968c9a 100644
--- a/src/callbacks_stage/subcell_bounds_check_2d.jl
+++ b/src/callbacks_stage/subcell_bounds_check_2d.jl
@@ -60,6 +60,20 @@
                 deviation_threaded[stride_size * Threads.threadid()] = deviation
             end
         end
+        for variable in limiter.positivity_variables_nonlinear
+            key = Symbol(string(variable), "_min")
+            deviation_threaded = idp_bounds_delta_local[key]
+            @threaded for element in eachelement(solver, cache)
+                deviation = deviation_threaded[stride_size * Threads.threadid()]
+                for j in eachnode(solver), i in eachnode(solver)
+                    var = variable(get_node_vars(u, equations, solver, i, j, element),
+                                   equations)
+                    deviation = max(deviation,
+                                    variable_bounds[key][i, j, element] - var)
+                end
+                deviation_threaded[stride_size * Threads.threadid()] = deviation
+            end
+        end
     end
 
     for (key, _) in idp_bounds_delta_local
@@ -92,6 +106,10 @@
                     print(f, ", ",
                           idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size])
                 end
+                for variable in limiter.positivity_variables_nonlinear
+                    print(f, ", ",
+                          idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size])
+                end
             end
             println(f)
         end
diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 3c6f759db2b..f5a632723cf 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -1632,6 +1632,18 @@ end
     return p
 end
 
+# Transformation from conservative variables u to d(p)/d(u)
+@inline function gradient_conservative(::typeof(pressure),
+                                       u, equations::CompressibleEulerEquations2D)
+    rho, rho_v1, rho_v2, rho_e = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v_square = v1^2 + v2^2
+
+    return (equations.gamma - 1.0) * SVector(0.5 * v_square, -v1, -v2, 1.0)
+end
+
 @inline function density_pressure(u, equations::CompressibleEulerEquations2D)
     rho, rho_v1, rho_v2, rho_e = u
     rho_times_p = (equations.gamma - 1) * (rho * rho_e - 0.5 * (rho_v1^2 + rho_v2^2))
@@ -1699,4 +1711,13 @@ end
 @inline function energy_internal(cons, equations::CompressibleEulerEquations2D)
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
+
+# State validation for Newton-bisection method of subcell IDP limiting
+@inline function Base.isvalid(u, equations::CompressibleEulerEquations2D)
+    p = pressure(u, equations)
+    if u[1] <= 0.0 || p <= 0.0
+        return false
+    end
+    return true
+end
 end # @muladd
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 7a3c326984d..c041bf117ba 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -376,6 +376,12 @@ of the correct length `nvariables(equations)`.
 """
 function energy_internal end
 
+# Default implementation of gradient for `variable`. Used for subcell limiting.
+# Implementing a gradient function for a specific variable improves the performance.
+@inline function gradient_conservative(variable, u, equations)
+    return ForwardDiff.gradient(x -> variable(x, equations), u)
+end
+
 ####################################################################################################
 # Include files with actual implementations for different systems of equations.
 
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index 43d1991e34b..4366cd32f08 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1118,6 +1118,20 @@ end
     return p
 end
 
+# Transformation from conservative variables u to d(p)/d(u)
+@inline function gradient_conservative(::typeof(pressure),
+                                       u, equations::IdealGlmMhdEquations2D)
+    rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v3 = rho_v3 / rho
+    v_square = v1^2 + v2^2 + v3^2
+
+    return (equations.gamma - 1.0) *
+           SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi)
+end
+
 @inline function density_pressure(u, equations::IdealGlmMhdEquations2D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
@@ -1384,6 +1398,15 @@ end
             cons[9]^2 / 2)
 end
 
+# State validation for Newton-bisection method of subcell IDP limiting
+@inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D)
+    p = pressure(u, equations)
+    if u[1] <= 0.0 || p <= 0.0
+        return false
+    end
+    return true
+end
+
 # Calculate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons'
 @inline function cross_helicity(cons, ::IdealGlmMhdEquations2D)
     return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1]
diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl
index 055e7ce24a4..e433c953779 100644
--- a/src/solvers/dgsem_tree/subcell_limiters.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters.jl
@@ -16,18 +16,28 @@ end
     SubcellLimiterIDP(equations::AbstractEquations, basis;
                       local_minmax_variables_cons = String[],
                       positivity_variables_cons = String[],
-                      positivity_correction_factor = 0.1)
+                      positivity_variables_nonlinear = [],
+                      positivity_correction_factor = 0.1,
+                      max_iterations_newton = 10,
+                      newton_tolerances = (1.0e-12, 1.0e-14),
+                      gamma_constant_newton = 2 * ndims(equations))
 
 Subcell invariant domain preserving (IDP) limiting used with [`VolumeIntegralSubcellLimiting`](@ref)
 including:
 - Local maximum/minimum Zalesak-type limiting for conservative variables (`local_minmax_variables_cons`)
-- Positivity limiting for conservative variables (`positivity_variables_cons`)
+- Positivity limiting for conservative variables (`positivity_variables_cons`) and nonlinear variables
+(`positivity_variables_nonlinear`)
 
 Conservative variables to be limited are passed as a vector of strings, e.g. `local_minmax_variables_cons = ["rho"]`
-and `positivity_variables_cons = ["rho"]`.
+and `positivity_variables_cons = ["rho"]`. For nonlinear variables the specific functions are
+passed in a vector, e.g. `positivity_variables_nonlinear = [pressure]`.
 
 The bounds are calculated using the low-order FV solution. The positivity limiter uses
 `positivity_correction_factor` such that `u^new >= positivity_correction_factor * u^FV`.
+The limiting of nonlinear variables uses a Newton-bisection method with a maximum of
+`max_iterations_newton` iterations, relative and absolute tolerances of `newton_tolerances`
+and a provisional update constant `gamma_constant_newton` (`gamma_constant_newton>=2*d`,
+where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) of Rueda-Ramírez et al. (2022).
 
 !!! note
     This limiter and the correction callback [`SubcellLimiterIDPCorrection`](@ref) only work together.
@@ -45,22 +55,32 @@ The bounds are calculated using the low-order FV solution. The positivity limite
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
 """
-struct SubcellLimiterIDP{RealT <: Real, Cache} <: AbstractSubcellLimiter
+struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear, Cache} <:
+       AbstractSubcellLimiter
     local_minmax::Bool
     local_minmax_variables_cons::Vector{Int}                   # Local mininum/maximum principles for conservative variables
     positivity::Bool
     positivity_variables_cons::Vector{Int}                     # Positivity for conservative variables
+    positivity_variables_nonlinear::LimitingVariablesNonlinear # Positivity for nonlinear variables
     positivity_correction_factor::RealT
     cache::Cache
+    max_iterations_newton::Int
+    newton_tolerances::Tuple{RealT, RealT}  # Relative and absolute tolerances for Newton's method
+    gamma_constant_newton::RealT            # Constant for the subcell limiting of convex (nonlinear) constraints
 end
 
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
 function SubcellLimiterIDP(equations::AbstractEquations, basis;
                            local_minmax_variables_cons = String[],
                            positivity_variables_cons = String[],
-                           positivity_correction_factor = 0.1)
+                           positivity_variables_nonlinear = [],
+                           positivity_correction_factor = 0.1,
+                           max_iterations_newton = 10,
+                           newton_tolerances = (1.0e-12, 1.0e-14),
+                           gamma_constant_newton = 2 * ndims(equations))
     local_minmax = (length(local_minmax_variables_cons) > 0)
-    positivity = (length(positivity_variables_cons) > 0)
+    positivity = (length(positivity_variables_cons) +
+                  length(positivity_variables_nonlinear) > 0)
 
     local_minmax_variables_cons_ = get_variable_index.(local_minmax_variables_cons,
                                                        equations)
@@ -80,13 +100,20 @@ function SubcellLimiterIDP(equations::AbstractEquations, basis;
             bound_keys = (bound_keys..., Symbol(string(v), "_min"))
         end
     end
+    for variable in positivity_variables_nonlinear
+        bound_keys = (bound_keys..., Symbol(string(variable), "_min"))
+    end
 
     cache = create_cache(SubcellLimiterIDP, equations, basis, bound_keys)
 
     SubcellLimiterIDP{typeof(positivity_correction_factor),
+                      typeof(positivity_variables_nonlinear),
                       typeof(cache)}(local_minmax, local_minmax_variables_cons_,
                                      positivity, positivity_variables_cons_,
-                                     positivity_correction_factor, cache)
+                                     positivity_variables_nonlinear,
+                                     positivity_correction_factor, cache,
+                                     max_iterations_newton, newton_tolerances,
+                                     gamma_constant_newton)
 end
 
 function Base.show(io::IO, limiter::SubcellLimiterIDP)
@@ -97,10 +124,15 @@ function Base.show(io::IO, limiter::SubcellLimiterIDP)
     if !(local_minmax || positivity)
         print(io, "No limiter selected => pure DG method")
     else
-        print(io, "limiter=(")
-        local_minmax && print(io, "min/max limiting, ")
-        positivity && print(io, "positivity")
-        print(io, "), ")
+        features = String[]
+        if local_minmax
+            push!(features, "local min/max")
+        end
+        if positivity
+            push!(features, "positivity")
+        end
+        join(io, features, ", ")
+        print(io, "Limiter=($features), ")
     end
     print(io, "Local bounds with FV solution")
     print(io, ")")
@@ -120,15 +152,15 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
             if local_minmax
                 setup = [
                     setup...,
-                    "" => "local maximum/minimum bounds for conservative variables $(limiter.local_minmax_variables_cons)",
+                    "" => "Local maximum/minimum limiting for conservative variables $(limiter.local_minmax_variables_cons)",
                 ]
             end
             if positivity
-                string = "positivity for conservative variables $(limiter.positivity_variables_cons)"
+                string = "Positivity limiting for conservative variables $(limiter.positivity_variables_cons) and $(limiter.positivity_variables_nonlinear)"
                 setup = [setup..., "" => string]
                 setup = [
                     setup...,
-                    "" => "   positivity correction factor = $(limiter.positivity_correction_factor)",
+                    "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)",
                 ]
             end
             setup = [
diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index 3d272359fe4..3f7954c8958 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -5,6 +5,10 @@
 @muladd begin
 #! format: noindent
 
+###############################################################################
+# IDP Limiting
+###############################################################################
+
 # this method is used when the limiter is constructed as for shock-capturing volume integrals
 function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquations{2},
                       basis::LobattoLegendreBasis, bound_keys)
@@ -66,6 +70,9 @@ function (limiter::SubcellLimiterIDP)(u::AbstractArray{<:Any, 4}, semi, dg::DGSE
     return nothing
 end
 
+###############################################################################
+# Calculation of local bounds using low-order FV solution
+
 @inline function calc_bounds_twosided!(var_min, var_max, variable, u, t, semi)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
     # Calc bounds inside elements
@@ -164,6 +171,9 @@ end
     return nothing
 end
 
+###############################################################################
+# Local minimum/maximum limiting
+
 @inline function idp_local_minmax!(alpha, limiter, u, t, dt, semi)
     for variable in limiter.local_minmax_variables_cons
         idp_local_minmax!(alpha, limiter, u, t, dt, semi, variable)
@@ -233,16 +243,36 @@ end
     return nothing
 end
 
+###############################################################################
+# Global positivity limiting
+
 @inline function idp_positivity!(alpha, limiter, u, dt, semi)
     # Conservative variables
     for variable in limiter.positivity_variables_cons
-        idp_positivity!(alpha, limiter, u, dt, semi, variable)
+        @trixi_timeit timer() "conservative variables" idp_positivity_conservative!(alpha,
+                                                                                    limiter,
+                                                                                    u,
+                                                                                    dt,
+                                                                                    semi,
+                                                                                    variable)
+    end
+
+    # Nonlinear variables
+    for variable in limiter.positivity_variables_nonlinear
+        @trixi_timeit timer() "nonlinear variables" idp_positivity_nonlinear!(alpha,
+                                                                              limiter,
+                                                                              u, dt,
+                                                                              semi,
+                                                                              variable)
     end
 
     return nothing
 end
 
-@inline function idp_positivity!(alpha, limiter, u, dt, semi, variable)
+###############################################################################
+# Global positivity limiting of conservative variables
+
+@inline function idp_positivity_conservative!(alpha, limiter, u, dt, semi, variable)
     mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
     (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
     (; inverse_weights) = dg.basis
@@ -256,7 +286,7 @@ end
         for j in eachnode(dg), i in eachnode(dg)
             var = u[variable, i, j, element]
             if var < 0
-                error("Safe $variable is not safe. element=$element, node: $i $j, value=$var")
+                error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.")
             end
 
             # Compute bound
@@ -302,4 +332,183 @@ end
 
     return nothing
 end
+
+@inline function idp_positivity_nonlinear!(alpha, limiter, u, dt, semi, variable)
+    _, equations, dg, cache = mesh_equations_solver_cache(semi)
+    (; positivity_correction_factor) = limiter
+
+    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
+    var_min = variable_bounds[Symbol(string(variable), "_min")]
+
+    @threaded for element in eachelement(dg, semi.cache)
+        inverse_jacobian = cache.elements.inverse_jacobian[element]
+        for j in eachnode(dg), i in eachnode(dg)
+            # Compute bound
+            u_local = get_node_vars(u, equations, dg, i, j, element)
+            var = variable(u_local, equations)
+            if var < 0
+                error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.")
+            end
+            var_min[i, j, element] = positivity_correction_factor * var
+
+            # Perform Newton's bisection method to find new alpha
+            newton_loops_alpha!(alpha, var_min[i, j, element], u_local, i, j, element,
+                                variable, initial_check_nonnegative_newton_idp,
+                                final_check_nonnegative_newton_idp, inverse_jacobian,
+                                dt, equations, dg, cache, limiter)
+        end
+    end
+
+    return nothing
+end
+
+@inline function newton_loops_alpha!(alpha, bound, u, i, j, element, variable,
+                                     initial_check, final_check, inverse_jacobian, dt,
+                                     equations, dg, cache, limiter)
+    (; inverse_weights) = dg.basis
+    (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
+
+    (; gamma_constant_newton) = limiter
+
+    # negative xi direction
+    antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[i] *
+                         get_node_vars(antidiffusive_flux1_R, equations, dg, i, j,
+                                       element)
+    newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check,
+                 equations, dt, limiter, antidiffusive_flux)
+
+    # positive xi direction
+    antidiffusive_flux = -gamma_constant_newton * inverse_jacobian *
+                         inverse_weights[i] *
+                         get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j,
+                                       element)
+    newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check,
+                 equations, dt, limiter, antidiffusive_flux)
+
+    # negative eta direction
+    antidiffusive_flux = gamma_constant_newton * inverse_jacobian * inverse_weights[j] *
+                         get_node_vars(antidiffusive_flux2_R, equations, dg, i, j,
+                                       element)
+    newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check,
+                 equations, dt, limiter, antidiffusive_flux)
+
+    # positive eta direction
+    antidiffusive_flux = -gamma_constant_newton * inverse_jacobian *
+                         inverse_weights[j] *
+                         get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1,
+                                       element)
+    newton_loop!(alpha, bound, u, i, j, element, variable, initial_check, final_check,
+                 equations, dt, limiter, antidiffusive_flux)
+
+    return nothing
+end
+
+@inline function newton_loop!(alpha, bound, u, i, j, element, variable, initial_check,
+                              final_check, equations, dt, limiter, antidiffusive_flux)
+    newton_reltol, newton_abstol = limiter.newton_tolerances
+
+    beta = 1 - alpha[i, j, element]
+
+    beta_L = 0 # alpha = 1
+    beta_R = beta # No higher beta (lower alpha) than the current one
+
+    u_curr = u + beta * dt * antidiffusive_flux
+
+    # If state is valid, perform initial check and return if correction is not needed
+    if isvalid(u_curr, equations)
+        goal = goal_function_newton_idp(variable, bound, u_curr, equations)
+
+        initial_check(bound, goal, newton_abstol) && return nothing
+    end
+
+    # Newton iterations
+    for iter in 1:(limiter.max_iterations_newton)
+        beta_old = beta
+
+        # If the state is valid, evaluate d(goal)/d(beta)
+        if isvalid(u_curr, equations)
+            dgoal_dbeta = dgoal_function_newton_idp(variable, u_curr, dt,
+                                                    antidiffusive_flux, equations)
+        else # Otherwise, perform a bisection step
+            dgoal_dbeta = 0
+        end
+
+        if dgoal_dbeta != 0
+            # Update beta with Newton's method
+            beta = beta - goal / dgoal_dbeta
+        end
+
+        # Check bounds
+        if (beta < beta_L) || (beta > beta_R) || (dgoal_dbeta == 0) || isnan(beta)
+            # Out of bounds, do a bisection step
+            beta = 0.5 * (beta_L + beta_R)
+            # Get new u
+            u_curr = u + beta * dt * antidiffusive_flux
+
+            # If the state is invalid, finish bisection step without checking tolerance and iterate further
+            if !isvalid(u_curr, equations)
+                beta_R = beta
+                continue
+            end
+
+            # Check new beta for condition and update bounds
+            goal = goal_function_newton_idp(variable, bound, u_curr, equations)
+            if initial_check(bound, goal, newton_abstol)
+                # New beta fulfills condition
+                beta_L = beta
+            else
+                # New beta does not fulfill condition
+                beta_R = beta
+            end
+        else
+            # Get new u
+            u_curr = u + beta * dt * antidiffusive_flux
+
+            # If the state is invalid, redefine right bound without checking tolerance and iterate further
+            if !isvalid(u_curr, equations)
+                beta_R = beta
+                continue
+            end
+
+            # Evaluate goal function
+            goal = goal_function_newton_idp(variable, bound, u_curr, equations)
+        end
+
+        # Check relative tolerance
+        if abs(beta_old - beta) <= newton_reltol
+            break
+        end
+
+        # Check absolute tolerance
+        if final_check(bound, goal, newton_abstol)
+            break
+        end
+    end
+
+    new_alpha = 1 - beta
+    if alpha[i, j, element] > new_alpha + newton_abstol
+        error("Alpha is getting smaller. old: $(alpha[i, j, element]), new: $new_alpha")
+    else
+        alpha[i, j, element] = new_alpha
+    end
+
+    return nothing
+end
+
+### Auxiliary routines for Newton's bisection method ###
+# Initial checks
+@inline initial_check_nonnegative_newton_idp(bound, goal, newton_abstol) = goal <= 0
+
+# Goal and d(Goal)d(u) function
+@inline goal_function_newton_idp(variable, bound, u, equations) = bound -
+                                                                  variable(u, equations)
+@inline function dgoal_function_newton_idp(variable, u, dt, antidiffusive_flux,
+                                           equations)
+    -dot(gradient_conservative(variable, u, equations), dt * antidiffusive_flux)
+end
+
+# Final checks
+@inline function final_check_nonnegative_newton_idp(bound, goal, newton_abstol)
+    (goal <= eps()) && (goal > -max(newton_abstol, abs(bound) * newton_abstol))
+end
 end # @muladd
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 61b5c54b5e9..b937abe92c0 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -581,6 +581,32 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"),
+                        l2=[
+                            0.42185634563805724,
+                            0.1686471269704017,
+                            0.18240674916968103,
+                            0.17858250604280654,
+                        ],
+                        linf=[
+                            1.7012978064377158,
+                            0.7149714986746726,
+                            0.5822547982757897,
+                            0.7300051017382696,
+                        ],
+                        tspan=(0.0, 2.0))
+    # 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)) < 15000
+    end
+end
+
 @trixi_testset "elixir_euler_colliding_flow.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"),
                         l2=[
diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl
index 953c077c0a3..1f8458075aa 100644
--- a/test/test_tree_2d_mhd.jl
+++ b/test/test_tree_2d_mhd.jl
@@ -332,24 +332,28 @@ end
 
 @trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"),
-                        l2=[2.9974425783503109e-02,
-                            7.2849646345685956e-02,
-                            7.2488477174662239e-02,
+                        l2=[
+                            3.2064026219236076e-02,
+                            7.2461094392606618e-02,
+                            7.2380202888062711e-02,
                             0.0000000000000000e+00,
-                            1.2507971380965512e+00,
-                            1.8929505145499678e-02,
-                            1.2218606317164420e-02,
+                            8.6293936673145932e-01,
+                            8.4091669534557805e-03,
+                            5.2156364913231732e-03,
                             0.0000000000000000e+00,
-                            3.0154796910479838e-03],
-                        linf=[3.2147382412340830e-01,
-                            1.3709471664007811e+00,
-                            1.3465154685288383e+00,
+                            2.0786952301129021e-04,
+                        ],
+                        linf=[
+                            3.8778760255775635e-01,
+                            9.4666683953698927e-01,
+                            9.4618924645661928e-01,
                             0.0000000000000000e+00,
-                            1.6051257523415284e+01,
-                            3.0564266749926644e-01,
-                            2.3908016329805595e-01,
+                            1.0980297261521951e+01,
+                            1.0264404591009069e-01,
+                            1.0655686942176350e-01,
                             0.0000000000000000e+00,
-                            1.3711262178549158e-01],
+                            6.1013422157115546e-03,
+                        ],
                         tspan=(0.0, 0.003))
     # Ensure that we do not have excessive memory allocations 
     # (e.g., from type instabilities) 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index e8a8effbe29..7943d952f71 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -416,7 +416,8 @@ end
     indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache")
     @test_nowarn show(stdout, indicator_hg)
 
-    limiter_idp = SubcellLimiterIDP(true, [1], true, [1], 0.1, "cache")
+    limiter_idp = SubcellLimiterIDP(true, [1], true, [1], ["variable"], 0.1, "cache", 1,
+                                    (1.0, 1.0), 1.0)
     @test_nowarn show(stdout, limiter_idp)
 
     # TODO: TrixiShallowWater: move unit test
@@ -1220,6 +1221,26 @@ end
     end
 end
 
+@testset "Consistency check for `gradient_conservative` routine" begin
+    # Set up conservative variables, equations
+    u = [
+        0.5011914484393387,
+        0.8829127712445113,
+        0.43024132987932817,
+        0.7560616633050348,
+    ]
+
+    equations = CompressibleEulerEquations2D(1.4)
+
+    # Define wrapper function for pressure in order to call default implementation
+    function pressure_test(u, equations)
+        return pressure(u, equations)
+    end
+
+    @test Trixi.gradient_conservative(pressure_test, u, equations) ≈
+          Trixi.gradient_conservative(pressure, u, equations)
+end
+
 @testset "Equivalent Fluxes" begin
     # Set up equations and dummy conservative variables state
     # Burgers' Equation