diff --git a/docs/src/index.md b/docs/src/index.md
index fe708585c..08359826a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,2 +1,9 @@
 # TurbulenceConvection.jl
 
+```@meta
+CurrentModule = TurbulenceConvection
+```
+
+```@docs
+TurbulenceConvection.ClimaParams
+```
diff --git a/integration_tests/utils/parameter_set.jl b/integration_tests/utils/parameter_set.jl
index ee6018782..5fb942524 100644
--- a/integration_tests/utils/parameter_set.jl
+++ b/integration_tests/utils/parameter_set.jl
@@ -22,6 +22,7 @@ function create_parameter_set(namelist)
         MSLP = 100000.0, # or grab from, e.g., namelist[""][...]
         τ_cond_evap = TC.parse_namelist(namelist, "microphysics", "τ_cond_evap"; default = 10.0),
         c_ε = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_factor"),
+        c_div = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "entrainment_massflux_div_factor"; default = 0.0),
         α_b = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_buoy_coeff1"),
         α_a = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_adv_coeff"),
         α_d = TC.parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "pressure_normalmode_drag_coeff"),
diff --git a/src/ClimaParams.jl b/src/ClimaParams.jl
new file mode 100644
index 000000000..9de861f1c
--- /dev/null
+++ b/src/ClimaParams.jl
@@ -0,0 +1,21 @@
+"""
+    ClimaParams
+
+This is a module for staging experimental
+clima parameters-- that we are unsure will
+be used in the final model. Once they are
+no longer experimental, they should be moved
+to CLIMAParameters.jl.
+
+TODO: move non-experimental parameters to CLIMAParameters
+"""
+module ClimaParams
+
+import CLIMAParameters
+const CP = CLIMAParameters
+const APS = CP.AbstractEarthParameterSet
+
+""" divergence factor for bubble case (zero otherwise) """
+entrainment_massflux_div_factor(ps::APS) = ps.nt.c_div
+
+end
diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
index 6d61f3371..46165db25 100644
--- a/src/TurbulenceConvection.jl
+++ b/src/TurbulenceConvection.jl
@@ -48,6 +48,10 @@ function parse_namelist(namelist, keys...; default = nothing, valid_options = no
     return param
 end
 
+include("ClimaParams.jl")
+import .ClimaParams
+const ICP = ClimaParams # internal clima parameters
+
 include("python_primitives.jl")
 include("parameters.jl")
 include("Grid.jl")
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index 6a6be5013..95fcc98ae 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -937,7 +937,6 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
                     w_min,
                     self.sorting_power,
                     self.detrainment_factor,
-                    self.entrainment_Mdiv_factor,
                     self.entrainment_sigma,
                     self.entrainment_scale,
                     self.updraft_mixing_frac,
diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
index 2dd360d3d..ae38a99d4 100644
--- a/src/closures/entr_detr.jl
+++ b/src/closures/entr_detr.jl
@@ -6,7 +6,6 @@
         w_min,
         β,
         c_δ,
-        c_div,
         c_μ,
         c_μ0,
         χ_upd,
@@ -21,7 +20,6 @@ Cohen et al. (JAMES, 2020), given:
  - `w_min`: minimum veritcal velocity
  - `β`: sorting power for moist mixing
  - `c_δ`: detrainment factor
- - `c_div`: divergence factor for bubble case (zero otherwise)
  - `c_μ`: logisitc function scale
  - `c_μ0`: logisitc function timescale
  - `χ_upd`: updraft mixing fraction
@@ -29,9 +27,10 @@ Cohen et al. (JAMES, 2020), given:
  - `c_εt`: turbulent entrainment factor
  - `εδ_model`: a [`MoistureDeficitEntr`](@ref)
 """
-function entr_detr(param_set, w_min, β, c_δ, c_div, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model::MoistureDeficitEntr)
+function entr_detr(param_set, w_min, β, c_δ, c_μ, c_μ0, χ_upd, c_λ, c_εt, εδ_model::MoistureDeficitEntr)
 
     c_ε = CPEDMF.c_ε(param_set)
+    c_div = ICP.entrainment_massflux_div_factor(param_set)
 
     if (εδ_model.q_liq_up + εδ_model.q_liq_en) == 0.0
         c_δ = 0.0
diff --git a/src/types.jl b/src/types.jl
index fe227042f..1aaac547a 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -737,7 +737,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
     extrapolate_buoyancy::Bool
     surface_area::Float64
     max_area::Float64
-    entrainment_Mdiv_factor::Float64
     updraft_mixing_frac::Float64
     entrainment_sigma::Float64
     entrainment_smin_tke_coeff::Float64
@@ -843,7 +842,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
         max_area = namelist["turbulence"]["EDMF_PrognosticTKE"]["max_area"]
         # entrainment parameters
         entrainment_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_factor"]
-        entrainment_Mdiv_factor = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_massflux_div_factor"]
         updraft_mixing_frac = namelist["turbulence"]["EDMF_PrognosticTKE"]["updraft_mixing_frac"]
         entrainment_sigma = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_sigma"]
         entrainment_smin_tke_coeff = namelist["turbulence"]["EDMF_PrognosticTKE"]["entrainment_smin_tke_coeff"]
@@ -960,7 +958,6 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
             extrapolate_buoyancy,
             surface_area,
             max_area,
-            entrainment_Mdiv_factor,
             updraft_mixing_frac,
             entrainment_sigma,
             entrainment_smin_tke_coeff,