diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index 06edfee53ca..f1daf89de54 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -229,6 +229,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜz = Fields.coordinate_field(Y.c).z
+    ᶜdz = Fields.Δz_field(axes(Y.c))
     (; params) = p
     (; dt) = p
     (; ᶜΦ, ᶜρ_ref) = p.core
@@ -295,6 +296,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
         ts_prev_level = Fields.field_values(Fields.level(ᶜts, i - 1))
         p_prev_level = Fields.field_values(Fields.level(ᶜp, i - 1))
         z_prev_level = Fields.field_values(Fields.level(ᶜz, i - 1))
+        dz_prev_level = Fields.field_values(Fields.level(ᶜdz, i - 1))
 
         local_geometry_prev_level = Fields.field_values(
             Fields.level(Fields.local_geometry_field(Y.c), i - 1),
@@ -368,6 +370,18 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 p.atmos.edmfx_entr_model,
             )
 
+            if i > 2
+                @. entrʲ_prev_level = limit_entrainment(
+                    entrʲ_prev_level,
+                    draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+                    get_physical_w(
+                        u³ʲ_prev_halflevel,
+                        local_geometry_prev_halflevel,
+                    ),
+                    dz_prev_level,
+                )
+            end
+
             # TODO: use updraft top instead of scale height
             @. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure(
                 params,
@@ -497,6 +511,18 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 p.atmos.edmfx_detr_model,
             )
 
+            if i > 2
+                @. detrʲ_prev_level = limit_entrainment(
+                    detrʲ_prev_level,
+                    draft_area(ρaʲ_prev_level, ρʲ_prev_level),
+                    get_physical_w(
+                        u³ʲ_prev_halflevel,
+                        local_geometry_prev_halflevel,
+                    ),
+                    dz_prev_level,
+                )
+            end
+
             ρaʲu³ʲ_data = p.scratch.temp_data_level_2
             ρaʲu³ʲ_datah_tot = ρaʲu³ʲ_dataq_tot = p.scratch.temp_data_level_3
 
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index 1a1279e24ef..8b2e76457b7 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -403,3 +403,9 @@ function edmfx_entr_detr_tendency!(
     end
     return nothing
 end
+
+limit_entrainment(entr::FT, a, dt) where {FT} =
+    max(min(entr, (1 - a) / max(a, eps(FT)) / dt), 0)
+limit_entrainment(entr, a, w, dz) = max(min(entr, w / dz), 0)
+limit_detrainment(detr, a, dt) = max(min(detr, 1 / dt), 0)
+limit_detrainment(detr, a, w, dz) = max(min(detr, w / dz), 0)