diff --git a/driver/Cases.jl b/driver/Cases.jl
index 7c348412f4..7fc22de98d 100644
--- a/driver/Cases.jl
+++ b/driver/Cases.jl
@@ -565,7 +565,7 @@ ForcingBase(case::TRMM_LBA, param_set::APS; kwargs...) =
 function surface_ref_state(::TRMM_LBA, param_set::APS, namelist)
     molmass_ratio = CPP.molmass_ratio(param_set)
     Pg = 991.3 * 100  #Pressure at ground
-    Tg = 296.85   # surface values for reference state (RS) which outputs p rho0 alpha0
+    Tg = 296.85   # surface values for reference state (RS) which outputs p, ρ
     pvg = TD.saturation_vapor_pressure(param_set, Tg, TD.Liquid())
     qtg = (1 / molmass_ratio) * pvg / (Pg - pvg) #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
@@ -647,7 +647,7 @@ ForcingBase(case::ARM_SGP, param_set::APS; kwargs...) =
 
 function surface_ref_state(::ARM_SGP, param_set::APS, namelist)
     Pg = 970.0 * 100 #Pressure at ground
-    Tg = 299.0   # surface values for reference state (RS) which outputs p rho0 alpha0
+    Tg = 299.0   # surface values for reference state (RS) which outputs  p, ρ
     qtg = 15.2 / 1000 #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
@@ -731,7 +731,7 @@ ForcingBase(case::GATE_III, param_set::APS; kwargs...) = ForcingBase(get_forcing
 
 function surface_ref_state(::GATE_III, param_set::APS, namelist)
     Pg = 1013.0 * 100  #Pressure at ground
-    Tg = 299.184   # surface values for reference state (RS) which outputs p rho0 alpha0
+    Tg = 299.184   # surface values for reference state (RS) which outputs p, ρ
     qtg = 16.5 / 1000 #Total water mixing ratio at surface
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
@@ -1028,7 +1028,7 @@ end
 function surface_ref_state(::SP, param_set::APS, namelist)
     Pg = 1.0e5  #Pressure at ground
     Tg = 300.0  #Temperature at ground
-    qtg = 1.0e-4   #Total water mixing ratio at TC. if set to 0, alpha0, rho0, p are NaN.
+    qtg = 1.0e-4   #Total water mixing ratio at TC. if set to 0, p, ρ are NaN.
     return TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
 end
 
diff --git a/src/closures/buoyancy_gradients.jl b/src/closures/buoyancy_gradients.jl
index b262071b82..fba87bc2f8 100644
--- a/src/closures/buoyancy_gradients.jl
+++ b/src/closures/buoyancy_gradients.jl
@@ -22,7 +22,7 @@ function buoyancy_gradients(
     phase_part = TD.PhasePartition(0.0, 0.0, 0.0) # assuming R_d = R_m
     Π = TD.exner_given_pressure(param_set, bg_model.p, phase_part)
 
-    ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p) * Π
+    ∂b∂θv = g * (R_d * bg_model.ρ / bg_model.p) * Π
 
     if bg_model.en_cld_frac > 0.0
         ts_sat = thermo_state_pθq(param_set, bg_model.p, bg_model.θ_liq_ice_sat, bg_model.qt_sat)
diff --git a/src/types.jl b/src/types.jl
index 995112c4e4..fe0fc3b9da 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -123,8 +123,8 @@ Base.@kwdef struct EnvBuoyGrad{FT, EBC <: AbstractEnvBuoyGradClosure}
     p::FT
     "cloud fraction"
     en_cld_frac::FT
-    "specific volume"
-    alpha0::FT
+    "density"
+    ρ::FT
 end
 function EnvBuoyGrad(::EBG; t_sat::FT, bg_kwargs...) where {FT <: Real, EBG <: AbstractEnvBuoyGradClosure}
     return EnvBuoyGrad{FT, EBG}(; t_sat, bg_kwargs...)
diff --git a/src/update_aux.jl b/src/update_aux.jl
index 48e498dc66..3b54c9c750 100644
--- a/src/update_aux.jl
+++ b/src/update_aux.jl
@@ -393,7 +393,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                 ∂θl∂z_sat = ∂θl∂z[k],
                 p = p_c[k],
                 en_cld_frac = aux_en.cloud_fraction[k],
-                alpha0 = 1 / ρ_c[k],
+                ρ = ρ_c[k],
             )
             bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)
 
@@ -410,7 +410,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas
                 ∂θl∂z_sat = ∂θl∂z_sat[k],
                 p = p_c[k],
                 en_cld_frac = aux_en.cloud_fraction[k],
-                alpha0 = 1 / ρ_c[k],
+                ρ = ρ_c[k],
             )
             bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)
         else