diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl
index 79f9c2e9403..e909c245b0b 100644
--- a/src/callbacks_stage/positivity_zhang_shu.jl
+++ b/src/callbacks_stage/positivity_zhang_shu.jl
@@ -40,9 +40,9 @@ end
 
 # Version used by the AMR callback
 function (limiter!::PositivityPreservingLimiterZhangShu)(u, mesh, equations, solver,
-                                                         cache)
+                                                         cache, args...)
     limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables, mesh, equations,
-                       solver, cache)
+                       solver, cache, args...)
 end
 
 # Iterate over tuples in a type-stable way using "lispy tuple programming",
@@ -54,21 +54,22 @@ end
 # compile times can increase otherwise - but a handful of elements per tuple
 # is definitely fine.
 function limiter_zhang_shu!(u, thresholds::NTuple{N, <:Real}, variables::NTuple{N, Any},
-                            mesh, equations, solver, cache) where {N}
+                            mesh, equations, solver, cache, args...) where {N}
     threshold = first(thresholds)
     remaining_thresholds = Base.tail(thresholds)
     variable = first(variables)
     remaining_variables = Base.tail(variables)
 
-    limiter_zhang_shu!(u, threshold, variable, mesh, equations, solver, cache)
+    limiter_zhang_shu!(u, threshold, variable, mesh, equations, solver, cache,
+                       args...)
     limiter_zhang_shu!(u, remaining_thresholds, remaining_variables, mesh, equations,
-                       solver, cache)
+                       solver, cache, args...)
     return nothing
 end
 
 # terminate the type-stable iteration over tuples
 function limiter_zhang_shu!(u, thresholds::Tuple{}, variables::Tuple{},
-                            mesh, equations, solver, cache)
+                            mesh, equations, solver, cache, args...)
     nothing
 end
 
diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl
index 7797eb95b09..21bdf1acded 100644
--- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl
+++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl
@@ -42,4 +42,107 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+
+# Modified version of the limiter used in the refinement step of the AMR callback.
+# To ensure admissibility after the refinement step, we compute a joint
+# limiting coefficient for all children elements and then limit against the
+# admissible mean value of the parent element.
+# This strategy is described in Remark 3 of the paper:
+# - Arpit Babbar, Praveen Chandrashekar (2025)
+#   Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with
+#   error based time stepping for hyperbolic conservation laws
+#   [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622)
+function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1},
+                            equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int}, u_mean_refined_elements)
+    @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`."
+
+    @threaded for i in eachindex(element_ids_new)
+        # Get the mean value from the parent element
+        u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i)
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+
+        theta = one(eltype(u))
+
+        # Iterate over the children of the current element to determine a joint limiting coefficient `theta`
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            # determine minimum value
+            value_min = typemax(eltype(u))
+            for i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, new_element_id)
+                value_min = min(value_min, variable(u_node, equations))
+            end
+
+            # detect if limiting is necessary
+            value_min < threshold || continue
+
+            theta = min(theta, (value_mean - threshold) / (value_mean - value_min))
+        end
+
+        theta < 1 || continue
+
+        # Make sure to really reach the threshold and not only by machine precision
+        theta -= eps(typeof(theta))
+
+        # Iterate again over the children to apply joint shifting
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            for i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, new_element_id)
+                set_node_vars!(u,
+                               theta * u_node + (1 - theta) * u_mean,
+                               equations, dg, i, new_element_id)
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Modified version of the limiter used in the coarsening step of the AMR callback.
+# To ensure admissibility after the coarsening step, we apply the limiter to
+# the coarsened elements.
+function limiter_zhang_shu!(u, threshold::Real, variable,
+                            mesh::AbstractMesh{1}, equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int})
+    @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
+
+    # Apply limiter to coarsened elements
+    @threaded for element in element_ids_new
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, element))
+        total_volume = zero(eltype(u))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            u_mean += u_node * weights[i]
+        end
+        # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
+        u_mean = u_mean / 2^ndims(mesh)
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl
index 813dd65878b..dc0730afad1 100644
--- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl
+++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl
@@ -47,4 +47,110 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+
+# Modified version of the limiter used in the refinement step of the AMR callback.
+# To ensure admissibility after the refinement step, we compute a joint
+# limiting coefficient for all children elements and then limit against the
+# admissible mean value of the parent element.
+# This strategy is described in Remark 3 of the paper:
+# - Arpit Babbar, Praveen Chandrashekar (2025)
+#   Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with
+#   error based time stepping for hyperbolic conservation laws
+#   [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622)
+function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{2},
+                            equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int}, u_mean_refined_elements)
+    @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`."
+
+    @threaded for i in eachindex(element_ids_new)
+        # Get the mean value from the parent element
+        u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i)
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+
+        theta = one(eltype(u))
+
+        # Iterate over the children of the current element to determine a joint limiting coefficient `theta`
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            # determine minimum value
+            value_min = typemax(eltype(u))
+            for j in eachnode(dg), i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, j, new_element_id)
+                value_min = min(value_min, variable(u_node, equations))
+            end
+
+            # detect if limiting is necessary
+            value_min < threshold || continue
+
+            theta = min(theta, (value_mean - threshold) / (value_mean - value_min))
+        end
+
+        theta < 1 || continue
+
+        # Make sure to really reach the threshold and not only by machine precision
+        theta -= eps(typeof(theta))
+
+        # Iterate again over the children to apply joint shifting
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            for j in eachnode(dg), i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, j, new_element_id)
+                set_node_vars!(u,
+                               theta * u_node + (1 - theta) * u_mean,
+                               equations, dg, i, j, new_element_id)
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Modified version of the limiter used in the coarsening step of the AMR callback.
+# To ensure admissibility after the coarsening step, we apply the limiter to
+# the coarsened elements.
+function limiter_zhang_shu!(u, threshold::Real, variable,
+                            mesh::AbstractMesh{2}, equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int})
+    @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
+
+    # Apply limiter to coarsened elements
+    @threaded for element in element_ids_new
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
+        total_volume = zero(eltype(u))
+        for j in eachnode(dg), i in eachnode(dg)
+            volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
+                                                           i, j, element)))
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            u_mean += u_node * weights[i] * weights[j] * volume_jacobian
+            total_volume += weights[i] * weights[j] * volume_jacobian
+        end
+        # normalize with the total volume
+        u_mean = u_mean / total_volume
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, j, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl
index 156abf35b4c..f8ce7c5d173 100644
--- a/src/callbacks_stage/positivity_zhang_shu_dg3d.jl
+++ b/src/callbacks_stage/positivity_zhang_shu_dg3d.jl
@@ -47,4 +47,110 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+
+# Modified version of the limiter used in the refinement step of the AMR callback.
+# To ensure admissibility after the refinement step, we compute a joint
+# limiting coefficient for all children elements and then limit against the
+# admissible mean value of the parent element.
+# This strategy is described in Remark 3 of the paper:
+# - Arpit Babbar, Praveen Chandrashekar (2025)
+#   Lax-Wendroff flux reconstruction on adaptive curvilinear meshes with
+#   error based time stepping for hyperbolic conservation laws
+#   [doi: 10.1016/j.jcp.2024.113622](https://doi.org/10.1016/j.jcp.2024.113622)
+function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{3},
+                            equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int}, u_mean_refined_elements)
+    @assert length(element_ids_new)==size(u_mean_refined_elements, 2) "The length of `element_ids_new` must match the second dimension of `u_mean_refined_elements`."
+
+    @threaded for i in eachindex(element_ids_new)
+        # Get the mean value from the parent element
+        u_mean = get_node_vars(u_mean_refined_elements, equations, dg, i)
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+
+        theta = one(eltype(u))
+
+        # Iterate over the children of the current element to determine a joint limiting coefficient `theta`
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            # determine minimum value
+            value_min = typemax(eltype(u))
+            for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id)
+                value_min = min(value_min, variable(u_node, equations))
+            end
+
+            # detect if limiting is necessary
+            value_min < threshold || continue
+
+            theta = min(theta, (value_mean - threshold) / (value_mean - value_min))
+        end
+
+        theta < 1 || continue
+
+        # Make sure to really reach the threshold and not only by machine precision
+        theta -= eps(typeof(theta))
+
+        # Iterate again over the children to apply joint shifting
+        for new_element_id in element_ids_new[i]:(element_ids_new[i] + 2^ndims(mesh) - 1)
+            for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+                u_node = get_node_vars(u, equations, dg, i, j, k, new_element_id)
+                set_node_vars!(u,
+                               theta * u_node + (1 - theta) * u_mean,
+                               equations, dg, i, j, k, new_element_id)
+            end
+        end
+    end
+
+    return nothing
+end
+
+# Modified version of the limiter used in the coarsening step of the AMR callback.
+# To ensure admissibility after the coarsening step, we apply the limiter to
+# the coarsened elements.
+function limiter_zhang_shu!(u, threshold::Real, variable,
+                            mesh::AbstractMesh{3}, equations, dg::DGSEM, cache,
+                            element_ids_new::Vector{Int})
+    @unpack weights = dg.basis
+    @unpack inverse_jacobian = cache.elements
+
+    # Apply limiter to coarsened elements
+    @threaded for element in element_ids_new
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
+        total_volume = zero(eltype(u))
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
+                                                           i, j, k, element)))
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian
+            total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian
+        end
+        # normalize with the total volume
+        u_mean = u_mean / total_volume
+
+        # We compute the value directly with the mean values, as we assume that
+        # Jensen's inequality holds (e.g. pressure for compressible Euler equations).
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            set_node_vars!(u, theta * u_node + (1 - theta) * u_mean,
+                           equations, dg, i, j, k, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl
index 196a5814f98..4ad8fac53d6 100644
--- a/src/callbacks_step/amr.jl
+++ b/src/callbacks_step/amr.jl
@@ -852,6 +852,38 @@ function original2refined(original_cell_ids, refined_original_cells, mesh)
     return shifted_cell_ids[original_cell_ids]
 end
 
+# Auxiliary function to compute the new element ids for refined elements
+# Used when applying positivity limiter after refinement step
+# Saves every first id of the 2^ndims child elements
+# All child elements can be addressed with element_ids_new[i]:(element_ids_new[i] + 2^ndims - 1)
+function compute_new_ids_refined_elements(elements_to_refine, mesh)
+    element_ids_new = copy(elements_to_refine)
+    for i in eachindex(element_ids_new)
+        # Each refined element increases all ids of the following elements by 2^ndims - 1
+        for j in (i + 1):length(element_ids_new)
+            element_ids_new[j] += 2^ndims(mesh) - 1
+        end
+    end
+
+    return element_ids_new
+end
+
+# Auxiliary function to compute the new element ids for removed elements
+# Used when applying positivity limiter after coarsening step
+function compute_new_ids_removed_elements(elements_to_remove, mesh)
+    @assert length(elements_to_remove) % (2^ndims(mesh))==0 "The length of `elements_to_remove` must be a multiple of 2^ndims(mesh)."
+
+    element_ids_new = zeros(Int, div(length(elements_to_remove), 2^ndims(mesh)))
+    for i in eachindex(element_ids_new)
+        # New element id is the id of the first child
+        # Additionally, all following ids decrease by (2^ndims - 1) per coarsened element
+        element_ids_new[i] = elements_to_remove[2^ndims(mesh) * (i - 1) + 1] -
+                             (2^ndims(mesh) - 1) * (i - 1)
+    end
+
+    return element_ids_new
+end
+
 """
     ControllerThreeLevel(semi, indicator; base_level=1,
                                           med_level=base_level, med_threshold=0.0,
diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl
index 4f202523188..492b1efb233 100644
--- a/src/callbacks_step/amr_dg1d.jl
+++ b/src/callbacks_step/amr_dg1d.jl
@@ -23,6 +23,27 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
     GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Compute mean values for the elements to be refined
+        # Only if limiter was passed
+        if limiter! !== nothing
+            (; weights) = dg.basis
+            u_mean_refined_elements = Matrix{eltype(u_ode)}(undef,
+                                                            nvariables(equations),
+                                                            length(elements_to_refine))
+            for element in eachindex(elements_to_refine)
+                old_element_id = elements_to_refine[element]
+                # compute mean value
+                u_mean = zero(get_node_vars(old_u, equations, dg, 1, old_element_id))
+                for i in eachnode(dg)
+                    u_node = get_node_vars(old_u, equations, dg, i, element)
+                    u_mean += u_node * weights[i]
+                end
+                # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
+                u_mean = u_mean / 2^ndims(mesh)
+                set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element)
+            end
+        end
+
         # Get new list of leaf cells
         leaf_cell_ids = local_leaf_cells(mesh.tree)
 
@@ -75,7 +96,12 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after refinement
+        element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new,
+                                                  u_mean_refined_elements)
     end
 
     return nothing
@@ -226,7 +252,11 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after coarsening
+        element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new)
     end
 
     return nothing
diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl
index 188d8a0ef2c..dce3b99d840 100644
--- a/src/callbacks_step/amr_dg2d.jl
+++ b/src/callbacks_step/amr_dg2d.jl
@@ -104,6 +104,32 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
     GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Compute mean values for the elements to be refined
+        # Only if limiter was passed
+        if limiter! !== nothing
+            (; weights) = dg.basis
+            u_mean_refined_elements = Matrix{eltype(u_ode)}(undef,
+                                                            nvariables(equations),
+                                                            length(elements_to_refine))
+            for element in eachindex(elements_to_refine)
+                old_element_id = elements_to_refine[element]
+                # compute mean value
+                u_mean = zero(get_node_vars(old_u, equations, dg, 1, 1, old_element_id))
+                total_volume = zero(eltype(old_u))
+                for j in eachnode(dg), i in eachnode(dg)
+                    volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian,
+                                                                   mesh, i, j,
+                                                                   old_element_id)))
+                    u_node = get_node_vars(old_u, equations, dg, i, j, old_element_id)
+                    u_mean += u_node * weights[i] * weights[j] * volume_jacobian
+                    total_volume += weights[i] * weights[j] * volume_jacobian
+                end
+                # normalize with the total volume
+                u_mean = u_mean / total_volume
+                set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element)
+            end
+        end
+
         if mesh isa P4estMesh
             # Loop over all elements in old container and scale the old solution by the Jacobian
             # prior to projection
@@ -176,7 +202,12 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after refinement
+        element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new,
+                                                  u_mean_refined_elements)
     end
 
     return nothing
@@ -382,7 +413,11 @@ function coarsen!(u_ode::AbstractVector, adaptor,
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after coarsening
+        element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new)
     end
 
     return nothing
diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl
index f0eb5e65771..b09845fad3d 100644
--- a/src/callbacks_step/amr_dg3d.jl
+++ b/src/callbacks_step/amr_dg3d.jl
@@ -33,6 +33,36 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM
     GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Compute mean values for the elements to be refined
+        # Only if limiter was passed
+        if limiter! !== nothing
+            (; weights) = dg.basis
+            u_mean_refined_elements = Matrix{eltype(u_ode)}(undef,
+                                                            nvariables(equations),
+                                                            length(elements_to_refine))
+            for element in eachindex(elements_to_refine)
+                old_element_id = elements_to_refine[element]
+                # compute mean value
+                u_mean = zero(get_node_vars(old_u, equations, dg, 1, 1, 1,
+                                            old_element_id))
+                total_volume = zero(eltype(old_u))
+                for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+                    volume_jacobian = abs(inv(get_inverse_jacobian(old_inverse_jacobian,
+                                                                   mesh, i, j, k,
+                                                                   old_element_id)))
+                    u_node = get_node_vars(old_u, equations, dg, i, j, k,
+                                           old_element_id)
+                    u_mean += u_node * weights[i] * weights[j] * weights[k] *
+                              volume_jacobian
+                    total_volume += weights[i] * weights[j] * weights[k] *
+                                    volume_jacobian
+                end
+                # normalize with the total volume
+                u_mean = u_mean / total_volume
+                set_node_vars!(u_mean_refined_elements, u_mean, equations, dg, element)
+            end
+        end
+
         if mesh isa P4estMesh
             # Loop over all elements in old container and scale the old solution by the Jacobian
             # prior to projection
@@ -108,7 +138,12 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after refinement
+        element_ids_new = compute_new_ids_refined_elements(elements_to_refine, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new,
+                                                  u_mean_refined_elements)
     end
 
     return nothing
@@ -301,7 +336,11 @@ function coarsen!(u_ode::AbstractVector, adaptor,
 
     # Apply the positivity limiter to the solution
     if limiter! !== nothing
-        limiter!(u, mesh, equations, dg, cache)
+        # Precompute list with new element ids after coarsening
+        element_ids_new = compute_new_ids_removed_elements(elements_to_remove, mesh)
+
+        @trixi_timeit timer() "limiter!" limiter!(u, mesh, equations, dg, cache,
+                                                  element_ids_new)
     end
 
     return nothing