diff --git a/src/callbacks_stage/subcell_bounds_check.jl b/src/callbacks_stage/subcell_bounds_check.jl
index 81e7007f5a3..d62cf501054 100644
--- a/src/callbacks_stage/subcell_bounds_check.jl
+++ b/src/callbacks_stage/subcell_bounds_check.jl
@@ -118,7 +118,7 @@ end
     if positivity
         for v in limiter.positivity_variables_cons
             println("$(variables[v]):\n- positivity: ",
-                    idp_bounds_delta[Symbol("$(v)_min")])
+                    idp_bounds_delta[Symbol("$(v)_min")][2])
         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 79b0925c444..2df10d0a757 100644
--- a/src/callbacks_stage/subcell_bounds_check_2d.jl
+++ b/src/callbacks_stage/subcell_bounds_check_2d.jl
@@ -12,35 +12,36 @@
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta) = limiter.cache
 
-    if save_errors
-        open("$output_directory/deviations.txt", "a") do f
-            print(f, iter, ", ", time)
-        end
-    end
     if positivity
         for v in limiter.positivity_variables_cons
             key = Symbol("$(v)_min")
-            deviation_min = zero(eltype(u))
+            deviation = idp_bounds_delta[key]
             for element in eachelement(solver, cache), j in eachnode(solver),
                 i in eachnode(solver)
 
                 var = u[v, i, j, element]
-                deviation_min = max(deviation_min,
-                                    variable_bounds[key][i, j, element] - var)
-            end
-            idp_bounds_delta[key] = max(idp_bounds_delta[key], deviation_min)
-            if save_errors
-                deviation_min_ = deviation_min
-                open("$output_directory/deviations.txt", "a") do f
-                    print(f, ", ", deviation_min_)
-                end
+                deviation[1] = max(deviation[1],
+                                   variable_bounds[key][i, j, element] - var)
             end
+            deviation[2] = max(deviation[2], deviation[1])
         end
     end
     if save_errors
+        # Print to output file
         open("$output_directory/deviations.txt", "a") do f
+            print(f, iter, ", ", time)
+            if positivity
+                for v in limiter.positivity_variables_cons
+                    key = Symbol("$(v)_min")
+                    print(f, ", ", idp_bounds_delta[key][1])
+                end
+            end
             println(f)
         end
+        # Reset first entries of idp_bounds_delta
+        for (key, _) in idp_bounds_delta
+            idp_bounds_delta[key][1] = zero(eltype(idp_bounds_delta[key][1]))
+        end
     end
 
     return nothing
diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index 690646e9b5a..4d32272898b 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -13,9 +13,9 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat
                                                                         nnodes(basis),
                                                                         bound_keys)
 
-    idp_bounds_delta = Dict{Symbol, real(basis)}()
+    idp_bounds_delta = Dict{Symbol, Vector{real(basis)}}()
     for key in bound_keys
-        idp_bounds_delta[key] = zero(real(basis))
+        idp_bounds_delta[key] = zeros(real(basis), 2)
     end
 
     return (; subcell_limiter_coefficients, idp_bounds_delta)