diff --git a/config/model_configs/prognostic_edmfx_bomex_box.yml b/config/model_configs/prognostic_edmfx_bomex_box.yml
index 2bf92df0ca..8d670d5c21 100644
--- a/config/model_configs/prognostic_edmfx_bomex_box.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_box.yml
@@ -23,7 +23,7 @@ y_elem: 2
 z_elem: 60
 z_stretch: false
 perturb_initstate: false
-dt: "20secs"
+dt: "10secs"
 t_end: "6hours"
 dt_save_state_to_disk: "30mins"
 toml: [toml/prognostic_edmfx_bomex_box.toml]
diff --git a/docs/Project.toml b/docs/Project.toml
index da563737fa..232eb5d8b8 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,8 +1,8 @@
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
-ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaAtmos = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
diff --git a/perf/Project.toml b/perf/Project.toml
index f4f3a45b7e..b3e8aa8ac6 100644
--- a/perf/Project.toml
+++ b/perf/Project.toml
@@ -5,7 +5,6 @@ ArtifactWrappers = "a14bc488-3040-4b00-9dc1-f6467924858a"
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaAnalysis = "29b5916a-a76c-4e73-9657-3c8fd22e65e6"
@@ -17,6 +16,7 @@ ClimaCorePlots = "cf7c7e5a-b407-4c48-9047-11a94a308626"
 ClimaCoreSpectra = "c2caaa1d-32ae-4754-ba0d-80e7561362e9"
 ClimaCoreTempestRemap = "d934ef94-cdd4-4710-83d6-720549644b70"
 ClimaCoreVTK = "c8b6d40d-e815-466f-95ae-c48aefa668fa"
+ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 ClimaTimeSteppers = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
index 02b15a1fb2..4b289ead88 100644
--- a/src/cache/diagnostic_edmf_precomputed_quantities.jl
+++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -483,6 +483,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 p_prev_level,
                 ρ_prev_level,
                 buoyancy_flux_sfc_halflevel,
+                ρaʲ_prev_level,
                 draft_area(ρaʲ_prev_level, ρʲ_prev_level),
                 get_physical_w(
                     u³ʲ_prev_halflevel,
@@ -498,6 +499,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 FT(0),
                 entrʲ_prev_level,
                 vert_div_level,
+                FT(0), # mass flux divergence is not implemented for diagnostic edmf
                 p.atmos.edmfx_detr_model,
             )
 
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index 51c2c92b2d..df85dbcb16 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -64,8 +64,6 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
     (; ᶜspecific, ᶜp, ᶜh_tot, ᶜK) = p.precomputed
     (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶠKᵥʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
-    ᶜinterp_lb = Operators.LeftBiasedF2C()
-    ᶜinterp_rb = Operators.RightBiasedF2C()
 
     for j in 1:n
         ᶜuʲ = ᶜuʲs.:($j)
@@ -195,6 +193,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     ᶜlg = Fields.local_geometry_field(Y.c)
 
     ᶜvert_div = p.scratch.ᶜtemp_scalar
+    ᶜmassflux_vert_div = p.scratch.ᶜtemp_scalar_2
     for j in 1:n
         # entrainment/detrainment
         @. ᶜentrʲs.:($$j) = entrainment(
@@ -219,6 +218,8 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             dt,
         )
         @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
+        @. ᶜmassflux_vert_div =
+            ᶜdivᵥ(ᶠinterp(Y.c.sgsʲs.:($$j).ρa) * ᶠu³ʲs.:($$j))
         @. ᶜdetrʲs.:($$j) = detrainment(
             params,
             ᶜz,
@@ -226,6 +227,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             ᶜp,
             Y.c.ρ,
             buoyancy_flux,
+            Y.c.sgsʲs.:($$j).ρa,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
             TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
@@ -235,6 +237,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             FT(0),
             ᶜentrʲs.:($$j),
             ᶜvert_div,
+            ᶜmassflux_vert_div,
             p.atmos.edmfx_detr_model,
         )
         @. ᶜdetrʲs.:($$j) = limit_detrainment(
diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl
index 7529890408..7fe01793d1 100644
--- a/src/parameters/Parameters.jl
+++ b/src/parameters/Parameters.jl
@@ -37,6 +37,7 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     detr_coeff::FT
     detr_buoy_coeff::FT
     detr_vertdiv_coeff::FT
+    detr_massflux_vertdiv_coeff::FT
     min_area_limiter_scale::FT
     min_area_limiter_power::FT
     max_area_limiter_scale::FT
diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
index d20e631567..f466cf1c74 100644
--- a/src/parameters/create_parameters.jl
+++ b/src/parameters/create_parameters.jl
@@ -24,6 +24,7 @@ function TurbulenceConvectionParameters(toml_dict::CP.AbstractTOMLDict)
         :mixing_length_smin_ub => :smin_ub,
         :EDMF_min_area => :min_area,
         :detr_vertdiv_coeff => :detr_vertdiv_coeff,
+        :detr_massflux_vertdiv_coeff => :detr_massflux_vertdiv_coeff,
         :max_area_limiter_power => :max_area_limiter_power,
         :min_area_limiter_power => :min_area_limiter_power,
         :pressure_normalmode_drag_coeff => :pressure_normalmode_drag_coeff,
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index d7863307d8..2cf7ccb3c8 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -114,38 +114,6 @@ function entrainment(
     return entr
 end
 
-function entrainment(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ::GeneralizedHarmonicsEntrainment,
-) where {FT}
-    turbconv_params = CAP.turbconv_params(params)
-    entr_inv_tau = CAP.entr_tau(turbconv_params)
-    entr_coeff = CAP.entr_coeff(turbconv_params)
-    min_area_limiter_scale = CAP.min_area_limiter_scale(turbconv_params)
-    min_area_limiter_power = CAP.min_area_limiter_power(turbconv_params)
-    a_min = CAP.min_area(turbconv_params)
-
-    min_area_limiter =
-        min_area_limiter_scale *
-        exp(-min_area_limiter_power * (max(ᶜaʲ, 0) - a_min))
-    entr =
-        entr_inv_tau + entr_coeff * abs(ᶜwʲ) / (ᶜz - z_sfc) + min_area_limiter
-
-    return entr * FT(2) * hm_limiter(ᶜaʲ)
-end
-
 """
    Return detrainment rate [1/s].
 
@@ -153,11 +121,12 @@ end
    - params set with model parameters
    - ᶜz, z_sfc, ᶜp, ᶜρ, - grid-scale height, surface height, grid-scale pressure and density
    - buoy_flux_surface - buoyancy flux at the surface
-   - ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft area, physical vertical velocity,
+   - ᶜρaʲ, ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft effective density, updraft area, physical vertical velocity,
                                    relative humidity and buoyancy
    - ᶜw⁰, ᶜRH⁰, ᶜbuoy⁰ - environment physical vertical velocity,
                               relative humidity and buoyancy
-   - dt - timestep
+   - ᶜentr - entrainment rate
+   - ᶜvert_div,ᶜmassflux_vert_div - vertical divergence, vertical mass flux divergence
 """
 function detrainment(
     params,
@@ -166,6 +135,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -175,6 +145,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -187,6 +158,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -196,6 +168,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
 
@@ -236,6 +209,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -245,6 +219,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
     turbconv_params = CAP.turbconv_params(params)
@@ -252,6 +227,8 @@ function detrainment(
     detr_coeff = CAP.detr_coeff(turbconv_params)
     detr_buoy_coeff = CAP.detr_buoy_coeff(turbconv_params)
     detr_vertdiv_coeff = CAP.detr_vertdiv_coeff(turbconv_params)
+    detr_massflux_vertdiv_coeff =
+        CAP.detr_massflux_vertdiv_coeff(turbconv_params)
     max_area_limiter_scale = CAP.max_area_limiter_scale(turbconv_params)
     max_area_limiter_power = CAP.max_area_limiter_power(turbconv_params)
     a_max = CAP.max_area(turbconv_params)
@@ -259,54 +236,20 @@ function detrainment(
     max_area_limiter =
         max_area_limiter_scale *
         exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    detr =
-        detr_inv_tau +
-        detr_coeff * abs(ᶜwʲ) +
-        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
-        max_area_limiter
-
-    return detr
-end
-
-function detrainment(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    ᶜentr::FT,
-    ᶜvert_div::FT,
-    ::GeneralizedHarmonicsDetrainment,
-) where {FT}
-    turbconv_params = CAP.turbconv_params(params)
-    detr_inv_tau = CAP.detr_tau(turbconv_params)
-    detr_coeff = CAP.detr_coeff(turbconv_params)
-    detr_buoy_coeff = CAP.detr_buoy_coeff(turbconv_params)
-    detr_vertdiv_coeff = CAP.detr_vertdiv_coeff(turbconv_params)
-    max_area_limiter_scale = CAP.max_area_limiter_scale(turbconv_params)
-    max_area_limiter_power = CAP.max_area_limiter_power(turbconv_params)
-    a_max = CAP.max_area(turbconv_params)
 
-    max_area_limiter =
-        max_area_limiter_scale *
-        exp(-max_area_limiter_power * (a_max - min(ᶜaʲ, 1)))
-    detr =
-        detr_inv_tau +
-        detr_coeff * abs(ᶜwʲ) +
-        detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
-        max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div +
-        max_area_limiter
+    if ᶜρaʲ <= 0
+        detr = 0
+    else
+        detr =
+            detr_inv_tau +
+            detr_coeff * abs(ᶜwʲ) +
+            detr_buoy_coeff * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
+            max(eps(FT), abs(ᶜwʲ - ᶜw⁰)) - detr_vertdiv_coeff * ᶜvert_div -
+            detr_massflux_vertdiv_coeff * min(ᶜmassflux_vert_div, 0) / ᶜρaʲ +
+            max_area_limiter
+    end
 
-    return detr * FT(2) * hm_limiter(ᶜaʲ)
+    return detr
 end
 
 function detrainment(
@@ -316,6 +259,7 @@ function detrainment(
     ᶜp::FT,
     ᶜρ::FT,
     buoy_flux_surface::FT,
+    ᶜρaʲ::FT,
     ᶜaʲ::FT,
     ᶜwʲ::FT,
     ᶜRHʲ::FT,
@@ -325,6 +269,7 @@ function detrainment(
     ᶜbuoy⁰::FT,
     ᶜentr::FT,
     ᶜvert_div::FT,
+    ᶜmassflux_vert_div::FT,
     ::ConstantAreaDetrainment,
 ) where {FT}
     detr = ᶜentr - ᶜvert_div
diff --git a/toml/prognostic_edmfx_bomex_box.toml b/toml/prognostic_edmfx_bomex_box.toml
index a3075317c4..7f9370356e 100644
--- a/toml/prognostic_edmfx_bomex_box.toml
+++ b/toml/prognostic_edmfx_bomex_box.toml
@@ -8,13 +8,13 @@ value = 1.0e-5
 value = 0.7
 
 [entr_inv_tau]
-value = 0
+value = 0.0002
 
 [entr_coeff]
-value = 0.1
+value = 0.2
 
 [min_area_limiter_scale]
-value = 0.001
+value = 0
 
 [min_area_limiter_power]
 value = 10
@@ -23,13 +23,16 @@ value = 10
 value = 0
 
 [detr_coeff]
-value = 2e-3
+value = 0
 
 [detr_buoy_coeff]
-value = 0.3
+value = 0
 
 [detr_vertdiv_coeff]
-value = 0.6
+value = 0
+
+[detr_massflux_vertdiv_coeff]
+value = 1
 
 [max_area_limiter_scale]
-value = 0.01
+value = 0
diff --git a/toml/prognostic_edmfx_dycoms_rf01_box.toml b/toml/prognostic_edmfx_dycoms_rf01_box.toml
index 47bd644a4c..b6b1616bbd 100644
--- a/toml/prognostic_edmfx_dycoms_rf01_box.toml
+++ b/toml/prognostic_edmfx_dycoms_rf01_box.toml
@@ -29,6 +29,9 @@ value = 1e-3
 value = 0.3
 
 [detr_vertdiv_coeff]
+value = 0
+
+[detr_massflux_vertdiv_coeff]
 value = 0.5
 
 [max_area_limiter_scale]
diff --git a/toml/prognostic_edmfx_simpleplume.toml b/toml/prognostic_edmfx_simpleplume.toml
index 4e967ad430..69ce495db5 100644
--- a/toml/prognostic_edmfx_simpleplume.toml
+++ b/toml/prognostic_edmfx_simpleplume.toml
@@ -1,47 +1,38 @@
 [EDMF_surface_area]
 value = 0.1
-type = "float"
 
 [EDMF_min_area]
-value = 1.0e-5
-type = "float" 
+value = 1.0e-5 
 
 [EDMF_max_area]
 value = 0.7
-type = "float"
 
 [entr_inv_tau]
 value = 0.0001
-type = "float"
 
 [entr_coeff]
 value = 0
-type = "float"
 
 [min_area_limiter_scale]
 value = 0
-type = "float"
 
 [min_area_limiter_power]
 value = 10
-type = "float"
 
 [detr_inv_tau]
 value = 0
-type = "float"
 
 [detr_coeff]
 value = 0
-type = "float"
 
 [detr_buoy_coeff]
-value = 0.3
-type = "float"
+value = 0
 
 [detr_vertdiv_coeff]
 value = 0
-type = "float"
+
+[detr_massflux_vertdiv_coeff]
+value = 1
 
 [max_area_limiter_scale]
-value = 0.01
-type = "float"
+value = 0