diff --git a/docs/src/APIs/Common/Thermodynamics.md b/docs/src/APIs/Common/Thermodynamics.md
index bb66d994305..7b5425e6b18 100644
--- a/docs/src/APIs/Common/Thermodynamics.md
+++ b/docs/src/APIs/Common/Thermodynamics.md
@@ -63,6 +63,7 @@ liquid_fraction
 liquid_ice_pottemp
 liquid_ice_pottemp_sat
 liquid_specific_humidity
+mixing_ratios
 moist_static_energy
 q_vap_saturation
 q_vap_saturation_liquid
@@ -74,6 +75,7 @@ saturation_adjustment
 saturation_excess
 saturation_vapor_pressure
 soundspeed_air
+shum_to_mixing_ratio
 specific_enthalpy
 specific_volume
 supersaturation
diff --git a/src/Common/Thermodynamics/relations.jl b/src/Common/Thermodynamics/relations.jl
index 61cd5bfb2bb..2450a362057 100644
--- a/src/Common/Thermodynamics/relations.jl
+++ b/src/Common/Thermodynamics/relations.jl
@@ -45,6 +45,8 @@ export liquid_fraction, PhasePartition_equil
 export dry_pottemp
 export virtual_pottemp
 export exner
+export shum_to_mixing_ratio
+export mixing_ratios
 export liquid_ice_pottemp
 export liquid_ice_pottemp_sat
 export relative_humidity
@@ -2185,6 +2187,44 @@ exner(ts::ThermodynamicState) = exner(
     PhasePartition(ts),
 )
 
+"""
+    shum_to_mixing_ratio(q, q_tot)
+
+Mixing ratio, from specific humidity
+ - `q` specific humidity
+ - `q_tot` total specific humidity
+"""
+function shum_to_mixing_ratio(q::FT, q_tot::FT) where {FT <: Real}
+    return q / (1 - q_tot)
+end
+
+"""
+    mixing_ratios(q::PhasePartition)
+
+Mixing ratios
+ - `r.tot` total mixing ratio
+ - `r.liq` liquid mixing ratio
+ - `r.ice` ice mixing ratio
+given a phase partition, `q`.
+"""
+function mixing_ratios(q::PhasePartition{FT}) where {FT <: Real}
+    return PhasePartition(
+        shum_to_mixing_ratio(q.tot, q.tot),
+        shum_to_mixing_ratio(q.liq, q.tot),
+        shum_to_mixing_ratio(q.ice, q.tot),
+    )
+end
+
+"""
+    mixing_ratios(ts::ThermodynamicState)
+
+Mixing ratios stored, in a phase partition, for
+ - total specific humidity
+ - liquid specific humidity
+ - ice specific humidity
+"""
+mixing_ratios(ts::ThermodynamicState) = mixing_ratios(PhasePartition(ts))
+
 """
     relative_humidity(param_set, T, p, phase_type, q::PhasePartition)
 
diff --git a/test/Atmos/EDMF/closures/entr_detr.jl b/test/Atmos/EDMF/closures/entr_detr.jl
index 7e071ea23c5..bf1a9405cc7 100644
--- a/test/Atmos/EDMF/closures/entr_detr.jl
+++ b/test/Atmos/EDMF/closures/entr_detr.jl
@@ -68,7 +68,13 @@ function entr_detr(
     lim_amp = entr.lim_amp
     w_min = entr.w_min
     # precompute vars
-    w_up_i = fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)
+    w_up_i = fix_void_up(
+        m.turbconv,
+        up[i].ρa*ρ_inv,
+        up[i].ρaw / up[i].ρa,
+        gm.ρu[3],
+        gm.ρu[3],
+    )
     sqrt_tke = sqrt(max(en.ρatke, 0) * ρ_inv / env.a)
     # ensure far from zero
     Δw = filter_w(w_up_i - env.w, w_min)
diff --git a/test/Atmos/EDMF/closures/mixing_length.jl b/test/Atmos/EDMF/closures/mixing_length.jl
index 0cd89ea2528..14262fd8e2a 100644
--- a/test/Atmos/EDMF/closures/mixing_length.jl
+++ b/test/Atmos/EDMF/closures/mixing_length.jl
@@ -84,7 +84,7 @@ function mixing_length(
     b = FT(0)
     a_up = vuntuple(i -> up[i].ρa * ρinv, N_up)
     w_up = vuntuple(N_up) do i
-        fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)
+        fix_void_up(m.turbconv, up[i].ρa/gm.ρ, up[i].ρaw / up[i].ρa)
     end
     b = sum(
         ntuple(N_up) do i
diff --git a/test/Atmos/EDMF/closures/pressure.jl b/test/Atmos/EDMF/closures/pressure.jl
index 29047b4b4df..6ef18b0e936 100644
--- a/test/Atmos/EDMF/closures/pressure.jl
+++ b/test/Atmos/EDMF/closures/pressure.jl
@@ -41,7 +41,13 @@ function perturbation_pressure(
     up_aux = aux.turbconv.updraft
     up_dif = diffusive.turbconv.updraft
 
-    w_up_i = fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)
+    w_up_i = fix_void_up(
+        m.turbconv,
+        up[i].ρa/state.ρ,
+        up[i].ρaw / up[i].ρa,
+        state.ρu[3],
+        state.ρu[3],
+    )
 
     nh_press_buoy = press.α_b * buoy.up[i]
     nh_pressure_adv = -press.α_a * w_up_i * up_dif[i].∇w[3]
diff --git a/test/Atmos/EDMF/edmf_kernels.jl b/test/Atmos/EDMF/edmf_kernels.jl
index c074ce42fe8..1f5b96a1594 100644
--- a/test/Atmos/EDMF/edmf_kernels.jl
+++ b/test/Atmos/EDMF/edmf_kernels.jl
@@ -301,6 +301,7 @@ function compute_gradient_argument!(
     gm = state
     up = state.turbconv.updraft
     en = state.turbconv.environment
+    ρ_inv = 1 / gm.ρ
 
     # Recover thermo states
     ts = recover_thermo_state_all(m, state, aux)
@@ -309,11 +310,16 @@ function compute_gradient_argument!(
     env = environment_vars(state, N_up)
 
     @unroll_map(N_up) do i
-        up_tf[i].w = fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)
+        up_tf[i].w = fix_void_up(
+            turbconv,
+            up[i].ρa * ρ_inv,
+            up[i].ρaw / up[i].ρa,
+            gm.ρu[3],
+            gm.ρu[3],
+        )
     end
     _grav::FT = grav(m.param_set)
 
-    ρ_inv = 1 / gm.ρ
     θ_liq_en = liquid_ice_pottemp(ts.en)
     q_tot_en = total_specific_humidity(ts.en)
 
@@ -413,14 +419,37 @@ end;
 
 function source(::EntrDetr{up_ρa{i}}, atmos, args) where {i}
     @unpack E_dyn, Δ_dyn, ρa_up = args.precomputed.turbconv
-    return fix_void_up(ρa_up[i], E_dyn[i] - Δ_dyn[i])
+    @unpack state = args
+    FT = eltype(args.state)
+    ρ_inv = 1 / args.state.ρ
+    return fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        E_dyn[i] - Δ_dyn[i],
+        FT(0),
+        FT(0),
+    )
 end
 
 function source(::EntrDetr{up_ρaw{i}}, atmos, args) where {i}
     @unpack E_dyn, Δ_dyn, E_trb, env, ρa_up, w_up = args.precomputed.turbconv
     up = args.state.turbconv.updraft
-    entr = fix_void_up(ρa_up[i], (E_dyn[i] + E_trb[i]) * env.w)
-    detr = fix_void_up(ρa_up[i], (Δ_dyn[i] + E_trb[i]) * w_up[i])
+    FT = eltype(args.state)
+    ρ_inv = 1 / args.state.ρ
+    entr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (E_dyn[i] + E_trb[i]) * env.w,
+        FT(0),
+        FT(0),
+    )
+    detr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (Δ_dyn[i] + E_trb[i]) * w_up[i],
+        FT(0),
+        FT(0),
+    )
 
     return entr - detr
 end
@@ -428,10 +457,23 @@ end
 function source(::EntrDetr{up_ρaθ_liq{i}}, atmos, args) where {i}
     @unpack E_dyn, Δ_dyn, E_trb, env, ρa_up, ts_en = args.precomputed.turbconv
     up = args.state.turbconv.updraft
+    FT = eltype(args.state)
+    ρ_inv = 1 / args.state.ρ
     θ_liq_en = liquid_ice_pottemp(ts_en)
-    entr = fix_void_up(ρa_up[i], (E_dyn[i] + E_trb[i]) * θ_liq_en)
-    detr =
-        fix_void_up(ρa_up[i], (Δ_dyn[i] + E_trb[i]) * up[i].ρaθ_liq / ρa_up[i])
+    entr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (E_dyn[i] + E_trb[i]) * θ_liq_en,
+        FT(0),
+        FT(0),
+    )
+    detr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (Δ_dyn[i] + E_trb[i]) * up[i].ρaθ_liq / ρa_up[i],
+        FT(0),
+        FT(0),
+    )
 
     return entr - detr
 end
@@ -439,10 +481,23 @@ end
 function source(::EntrDetr{up_ρaq_tot{i}}, atmos, args) where {i}
     @unpack E_dyn, Δ_dyn, E_trb, env, ρa_up, ts_en = args.precomputed.turbconv
     up = args.state.turbconv.updraft
+    FT = eltype(args.state)
+    ρ_inv = 1 / args.state.ρ
     q_tot_en = total_specific_humidity(ts_en)
-    entr = fix_void_up(ρa_up[i], (E_dyn[i] + E_trb[i]) * q_tot_en)
-    detr =
-        fix_void_up(ρa_up[i], (Δ_dyn[i] + E_trb[i]) * up[i].ρaq_tot / ρa_up[i])
+    entr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (E_dyn[i] + E_trb[i]) * q_tot_en,
+        FT(0),
+        FT(0),
+    )
+    detr = fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        (Δ_dyn[i] + E_trb[i]) * up[i].ρaq_tot / ρa_up[i],
+        FT(0),
+        FT(0),
+    )
 
     return entr - detr
 end
@@ -453,16 +508,20 @@ function source(::EntrDetr{en_ρatke}, atmos, args)
     up = state.turbconv.updraft
     en = state.turbconv.environment
     gm = state
+    FT = eltype(args.state)
     N_up = n_updrafts(atmos.turbconv)
     ρ_inv = 1 / gm.ρ
     tke_en = enforce_positivity(en.ρatke) * ρ_inv / env.a
 
     entr_detr = vuntuple(N_up) do i
         fix_void_up(
-            ρa_up[i],
+            atmos.turbconv,
+            ρa_up[i] * ρ_inv,
             E_trb[i] * (env.w - gm.ρu[3] * ρ_inv) * (env.w - w_up[i]) -
             (E_dyn[i] + E_trb[i]) * tke_en +
             Δ_dyn[i] * (w_up[i] - env.w) * (w_up[i] - env.w) / 2,
+            FT(0),
+            FT(0),
         )
     end
     return sum(entr_detr)
@@ -471,16 +530,19 @@ end
 function source(::EntrDetr{en_ρaθ_liq_cv}, atmos, args)
     @unpack E_dyn, Δ_dyn, E_trb, ρa_up, ts_en = args.precomputed.turbconv
     @unpack state = args
+    FT = eltype(state)
     ts_gm = args.precomputed.ts
     up = state.turbconv.updraft
     en = state.turbconv.environment
     N_up = n_updrafts(atmos.turbconv)
     θ_liq = liquid_ice_pottemp(ts_gm)
     θ_liq_en = liquid_ice_pottemp(ts_en)
+    ρ_inv = 1 / args.state.ρ
 
     entr_detr = vuntuple(N_up) do i
         fix_void_up(
-            ρa_up[i],
+            atmos.turbconv,
+            ρa_up[i] * ρ_inv,
             Δ_dyn[i] *
             (up[i].ρaθ_liq / ρa_up[i] - θ_liq_en) *
             (up[i].ρaθ_liq / ρa_up[i] - θ_liq_en) +
@@ -491,6 +553,8 @@ function source(::EntrDetr{en_ρaθ_liq_cv}, atmos, args)
             (θ_liq_en - θ_liq) *
             (θ_liq_en - up[i].ρaθ_liq / ρa_up[i]) -
             (E_dyn[i] + E_trb[i]) * en.ρaθ_liq_cv,
+            FT(0),
+            FT(0),
         )
     end
     return sum(entr_detr)
@@ -510,7 +574,8 @@ function source(::EntrDetr{en_ρaq_tot_cv}, atmos, args)
 
     entr_detr = vuntuple(N_up) do i
         fix_void_up(
-            ρa_up[i],
+            atmos.turbconv,
+            ρa_up[i] * ρ_inv,
             Δ_dyn[i] *
             (up[i].ρaq_tot / ρa_up[i] - q_tot_en) *
             (up[i].ρaq_tot / ρa_up[i] - q_tot_en) +
@@ -521,6 +586,8 @@ function source(::EntrDetr{en_ρaq_tot_cv}, atmos, args)
             (q_tot_en - ρq_tot * ρ_inv) *
             (q_tot_en - up[i].ρaq_tot / ρa_up[i]) -
             (E_dyn[i] + E_trb[i]) * en.ρaq_tot_cv,
+            FT(0),
+            FT(0),
         )
     end
     return sum(entr_detr)
@@ -543,7 +610,8 @@ function source(::EntrDetr{en_ρaθ_liq_q_tot_cv}, atmos, args)
 
     entr_detr = vuntuple(N_up) do i
         fix_void_up(
-            ρa_up[i],
+            atmos.turbconv,
+            ρa_up[i] * ρ_inv,
             Δ_dyn[i] *
             (up[i].ρaθ_liq / ρa_up[i] - θ_liq_en) *
             (up[i].ρaq_tot / ρa_up[i] - q_tot_en) +
@@ -554,6 +622,8 @@ function source(::EntrDetr{en_ρaθ_liq_q_tot_cv}, atmos, args)
             (q_tot_en - ρq_tot * ρ_inv) *
             (θ_liq_en - up[i].ρaθ_liq / ρa_up[i]) -
             (E_dyn[i] + E_trb[i]) * en.ρaθ_liq_q_tot_cv,
+            FT(0),
+            FT(0),
         )
     end
     return sum(entr_detr)
@@ -562,9 +632,17 @@ end
 function source(::PressSource{en_ρatke}, atmos, args)
     @unpack env, ρa_up, dpdz, w_up = args.precomputed.turbconv
     up = args.state.turbconv.updraft
+    FT = eltype(args.state)
+    ρ_inv = 1 / args.state.ρ
     N_up = n_updrafts(atmos.turbconv)
     press_tke = vuntuple(N_up) do i
-        fix_void_up(ρa_up[i], ρa_up[i] * (w_up[i] - env.w) * dpdz[i])
+        fix_void_up(
+            atmos.turbconv,
+            ρa_up[i] * ρ_inv,
+            ρa_up[i] * (w_up[i] - env.w) * dpdz[i],
+            FT(0),
+            FT(0),
+        )
     end
     return sum(press_tke)
 end
@@ -688,7 +766,7 @@ function compute_ρa_up(atmos, state, aux)
     a_max = turbconv.subdomains.a_max
     # in future GCM implementations we need to think about grid mean advection
     ρa_up = vuntuple(N_up) do i
-        gm.ρ * enforce_unit_bounds(up[i].ρa / gm.ρ, a_min, a_max)
+        fix_void_up(turbconv, up[i].ρa / gm.ρ, up[i].ρa)
     end
     return ρa_up
 end
@@ -696,32 +774,66 @@ end
 function flux(::Advect{up_ρa{i}}, atmos, args) where {i}
     @unpack state, aux = args
     @unpack ρa_up = args.precomputed.turbconv
+    FT = eltype(state)
     up = state.turbconv.updraft
+    ρ_inv = 1 / args.state.ρ
     ẑ = vertical_unit_vector(atmos, aux)
-    return fix_void_up(ρa_up[i], up[i].ρaw) * ẑ
+    return fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        up[i].ρaw,
+        FT(0),
+        args.state.ρu[3],
+    ) * ẑ
 end
 function flux(::Advect{up_ρaw{i}}, atmos, args) where {i}
     @unpack state, aux = args
     @unpack ρa_up, w_up = args.precomputed.turbconv
+    FT = eltype(state)
     up = state.turbconv.updraft
+    ρ_inv = 1 / args.state.ρ
     ẑ = vertical_unit_vector(atmos, aux)
-    return fix_void_up(ρa_up[i], up[i].ρaw * w_up[i]) * ẑ
+    return fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        up[i].ρaw * w_up[i],
+        FT(0),
+        args.state.ρu[3] * args.state.ρu[3] / args.state.ρ,
+    ) * ẑ
 
 end
 function flux(::Advect{up_ρaθ_liq{i}}, atmos, args) where {i}
     @unpack state, aux = args
     @unpack ρa_up, w_up = args.precomputed.turbconv
+    ts_gm = args.precomputed.ts
+    FT = eltype(state)
     up = state.turbconv.updraft
+    ρ_inv = 1 / args.state.ρ
     ẑ = vertical_unit_vector(atmos, aux)
-    return fix_void_up(ρa_up[i], w_up[i] * up[i].ρaθ_liq) * ẑ
+    θ_liq = liquid_ice_pottemp(ts_gm)
+    return fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        up[i].ρaw * w_up[i],
+        FT(0),
+        args.state.ρu[3] * θ_liq,
+    ) * ẑ
 
 end
 function flux(::Advect{up_ρaq_tot{i}}, atmos, args) where {i}
     @unpack state, aux = args
     @unpack ρa_up, w_up = args.precomputed.turbconv
     up = state.turbconv.updraft
+    FT = eltype(state)
+    ρ_inv = 1 / args.state.ρ
     ẑ = vertical_unit_vector(atmos, aux)
-    return fix_void_up(ρa_up[i], w_up[i] * up[i].ρaq_tot) * ẑ
+    return fix_void_up(
+        atmos.turbconv,
+        ρa_up[i] * ρ_inv,
+        up[i].ρaq_tot * w_up[i],
+        FT(0),
+        args.state.ρu[3] * args.state.moisture.ρq_tot / args.state.ρ,
+    ) * ẑ
 
 end
 
@@ -796,7 +908,13 @@ function precompute(::EDMF, bl, args, ts, ::Flux{FirstOrder})
     N_up = n_updrafts(bl.turbconv)
     ρ_inv = 1 / gm.ρ
     w_up = vuntuple(N_up) do i
-        fix_void_up(ρa_up[i], up[i].ρaw / ρa_up[i])
+        fix_void_up(
+            bl.turbconv,
+            ρa_up[i] * ρ_inv,
+            up[i].ρaw / ρa_up[i],
+            state.ρu[3],
+            state.ρu[3],
+        )
     end
 
     θ_liq_up = vuntuple(N_up) do i
@@ -915,6 +1033,7 @@ function precompute(::EDMF, bl, args, ts, ::Source)
     ts_gm = ts
     gm = state
     up = state.turbconv.updraft
+    ρ_inv = 1 / args.state.ρ
     N_up = n_updrafts(bl.turbconv)
     # Get environment variables
     env = environment_vars(state, N_up)
@@ -940,7 +1059,13 @@ function precompute(::EDMF, bl, args, ts, ::Source)
     )
 
     w_up = vuntuple(N_up) do i
-        fix_void_up(ρa_up[i], up[i].ρaw / ρa_up[i])
+        fix_void_up(
+            bl.turbconv,
+            ρa_up[i] * ρ_inv,
+            up[i].ρaw / ρa_up[i],
+            state.ρu[3],
+            state.ρu[3],
+        )
     end
 
     en = state.turbconv.environment
@@ -991,10 +1116,15 @@ function flux(::SGSFlux{Energy}, atmos, args)
     massflux_e = sum(
         ntuple(N_up) do i
             fix_void_up(
-                ρa_up[i],
+                atmos.turbconv,
+                ρa_up[i] * ρ_inv,
+                atmos.turbconv,
+                ρa_up[i] * ρ_inv,
                 ρa_up[i] *
                 (gm.ρe * ρ_inv - e_tot_up[i]) *
                 (gm.ρu[3] * ρ_inv - ρaw_up[i] / ρa_up[i]),
+                fallback_low = FT(0),
+                fallback_high = FT(0),
             )
         end,
     )
@@ -1019,10 +1149,13 @@ function flux(::SGSFlux{TotalMoisture}, atmos, args)
     massflux_q_tot = sum(
         ntuple(N_up) do i
             fix_void_up(
-                ρa_up[i],
+                atmos.turbconv,
+                ρa_up[i] * ρ_inv,
                 ρa_up[i] *
                 (ρq_tot * ρ_inv - ρaq_tot_up[i] / ρa_up[i]) *
                 (ρu_gm_tup[3] * ρ_inv - ρaw_up[i] / ρa_up[i]),
+                fallback_low = FT(0),
+                fallback_high = FT(0),
             )
         end,
     )
@@ -1046,10 +1179,13 @@ function flux(::SGSFlux{Momentum}, atmos, args)
     massflux_w = sum(
         ntuple(N_up) do i
             fix_void_up(
-                ρa_up[i],
+                atmos.turbconv,
+                ρa_up[i] * ρ_inv,
                 ρa_up[i] *
                 (ρu_gm_tup[3] * ρ_inv - ρaw_up[i] / ρa_up[i]) *
                 (ρu_gm_tup[3] * ρ_inv - ρaw_up[i] / ρa_up[i]),
+                fallback_low = FT(0),
+                fallback_high = FT(0),
             )
         end,
     )
diff --git a/test/Atmos/EDMF/helper_funcs/nondimensional_exchange_functions.jl b/test/Atmos/EDMF/helper_funcs/nondimensional_exchange_functions.jl
index c12e0ac648c..777e44fc391 100644
--- a/test/Atmos/EDMF/helper_funcs/nondimensional_exchange_functions.jl
+++ b/test/Atmos/EDMF/helper_funcs/nondimensional_exchange_functions.jl
@@ -47,7 +47,13 @@ function nondimensional_exchange_functions(
     N_up = n_updrafts(m.turbconv)
     ρinv = 1 / gm.ρ
     a_up_i = up[i].ρa * ρinv
-    w_up_i = fix_void_up(up[i].ρa, up[i].ρaw / up[i].ρa)
+    w_up_i = fix_void_up(
+        m.turbconv,
+        up[i].ρa,
+        up[i].ρaw / up[i].ρa,
+        gm.ρu[3],
+        gm.ρu[3],
+    )
 
     # thermodynamic variables
     RH_up = relative_humidity(ts_up[i])
diff --git a/test/Atmos/EDMF/helper_funcs/subdomain_thermo_states.jl b/test/Atmos/EDMF/helper_funcs/subdomain_thermo_states.jl
index 0b722d308c2..7c573cc02fb 100644
--- a/test/Atmos/EDMF/helper_funcs/subdomain_thermo_states.jl
+++ b/test/Atmos/EDMF/helper_funcs/subdomain_thermo_states.jl
@@ -153,8 +153,13 @@ function new_thermo_state_up(
     ts_up = vuntuple(N_up) do i
         ρa_up = up[i].ρa
         ρaθ_liq_up = up[i].ρaθ_liq
-        θ_liq_up =
-            fix_void_up(ρa_up, ρaθ_liq_up / ρa_up, liquid_ice_pottemp(ts))
+        θ_liq_up = fix_void_up(
+            m.turbconv,
+            ρa_up/state.ρ,
+            ρaθ_liq_up / ρa_up,
+            liquid_ice_pottemp(ts),
+            liquid_ice_pottemp(ts),
+        )
         PhaseDry_pθ(m.param_set, p, θ_liq_up)
     end
     return ts_up
@@ -177,12 +182,19 @@ function new_thermo_state_up(
         ρa_up = up[i].ρa
         ρaθ_liq_up = up[i].ρaθ_liq
         ρaq_tot_up = up[i].ρaq_tot
-        θ_liq_up =
-            fix_void_up(ρa_up, ρaθ_liq_up / ρa_up, liquid_ice_pottemp(ts))
+        θ_liq_up = fix_void_up(
+            m.turbconv,
+            ρa_up/state.ρ,
+            ρaθ_liq_up / ρa_up,
+            liquid_ice_pottemp(ts),
+            liquid_ice_pottemp(ts),
+        )
         q_tot_up = fix_void_up(
-            ρa_up,
+            m.turbconv,
+            ρa_up/state.ρ,
             ρaq_tot_up / ρa_up,
             total_specific_humidity(ts),
+            total_specific_humidity(ts),
         )
         PhaseEquil_pθq(m.param_set, p, θ_liq_up, q_tot_up)
     end
diff --git a/test/Atmos/EDMF/helper_funcs/utility_funcs.jl b/test/Atmos/EDMF/helper_funcs/utility_funcs.jl
index bd5542b2d7f..d1f28b6af08 100644
--- a/test/Atmos/EDMF/helper_funcs/utility_funcs.jl
+++ b/test/Atmos/EDMF/helper_funcs/utility_funcs.jl
@@ -39,7 +39,27 @@ enforce_positivity(x::FT) where {FT} = max(x, FT(0))
 Substitute value by a consistent fallback in case of
 negligible area fraction (void updraft).
 """
-function fix_void_up(ρa_up_i::FT, val::FT, fallback = FT(0)) where {FT}
-    tol = sqrt(eps(FT))
-    return ρa_up_i > tol ? val : fallback
+# function fix_void_up(a_up_i::FT, val::FT, fallback = FT(0)) where {FT}
+#     tol = sqrt(eps(FT))
+#     return a_up_i > tol ? val : fallback
+# end
+
+function fix_void_up(
+    turbconv::EDMF,
+    a_up_i::FT,
+    val::FT,
+    fallback_low = FT(0),
+    fallback_high = FT(1),
+) where {FT}
+    # tol = sqrt(eps(FT))
+    a_min = turbconv.subdomains.a_min
+    a_max = turbconv.subdomains.a_max
+    if a_up_i > a_max
+        output = fallback_high
+    elseif a_up_i < a_min
+        output = fallback_low
+    else
+        output = val
+    end
+    return output
 end
diff --git a/test/Atmos/EDMF/report_mse_bomex.jl b/test/Atmos/EDMF/report_mse_bomex.jl
index 7524686efe8..13e73ec950d 100644
--- a/test/Atmos/EDMF/report_mse_bomex.jl
+++ b/test/Atmos/EDMF/report_mse_bomex.jl
@@ -16,13 +16,13 @@ best_mse[:Bomex]["ρ"] = 3.4917543567416755e-02
 best_mse[:Bomex]["ρu[1]"] = 3.0715061616086027e+03
 best_mse[:Bomex]["ρu[2]"] = 1.2895273328644972e-03
 best_mse[:Bomex]["moisture.ρq_tot"] = 4.1330591681441348e-02
-best_mse[:Bomex]["turbconv.environment.ρatke"] = 6.6415719930880925e+02
-best_mse[:Bomex]["turbconv.environment.ρaθ_liq_cv"] = 8.5667223192888514e+01
+best_mse[:Bomex]["turbconv.environment.ρatke"] = 6.6868199378725114e+02
+best_mse[:Bomex]["turbconv.environment.ρaθ_liq_cv"] = 8.5667224210019953e+01
 best_mse[:Bomex]["turbconv.environment.ρaq_tot_cv"] = 1.6435555167634794e+02
-best_mse[:Bomex]["turbconv.updraft[1].ρa"] = 7.9564915645201182e+01
-best_mse[:Bomex]["turbconv.updraft[1].ρaw"] = 8.4288782126742318e-02
-best_mse[:Bomex]["turbconv.updraft[1].ρaθ_liq"] = 9.0095910670762631e+00
-best_mse[:Bomex]["turbconv.updraft[1].ρaq_tot"] = 1.0768554319447651e+01
+best_mse[:Bomex]["turbconv.updraft[1].ρa"] = 7.9962915358154191e+01
+best_mse[:Bomex]["turbconv.updraft[1].ρaw"] = 8.5167510817010642e-02
+best_mse[:Bomex]["turbconv.updraft[1].ρaθ_liq"] = 9.0184987300428983e+00
+best_mse[:Bomex]["turbconv.updraft[1].ρaq_tot"] = 1.0779279608651940e+01
 #! format: on
 
 sufficient_mse(computed_mse, best_mse) = computed_mse <= best_mse + sqrt(eps())
diff --git a/test/Common/Thermodynamics/runtests.jl b/test/Common/Thermodynamics/runtests.jl
index 90c55842e17..b97766ce5fe 100644
--- a/test/Common/Thermodynamics/runtests.jl
+++ b/test/Common/Thermodynamics/runtests.jl
@@ -430,6 +430,20 @@ end
     q_tot = FT(0.23)
     @test TD.exner_given_pressure(param_set, p, PhasePartition(q_tot)) isa
           typeof(p)
+
+    q_tot = 0.1
+    q_liq = 0.05
+    q_ice = 0.01
+    mr = shum_to_mixing_ratio(q_tot, q_tot)
+    @test mr == q_tot / (1 - q_tot)
+    mr = shum_to_mixing_ratio(q_liq, q_tot)
+    @test mr == q_liq / (1 - q_tot)
+
+    q = PhasePartition(q_tot, q_liq, q_ice)
+    mrs = mixing_ratios(q)
+    @test mrs.tot == q_tot / (1 - q_tot)
+    @test mrs.liq == q_liq / (1 - q_tot)
+    @test mrs.ice == q_ice / (1 - q_tot)
 end
 
 
@@ -1087,6 +1101,10 @@ end
     @test PhasePartition(ts_eq).liq ≈ PhasePartition(ts_dry).liq
     @test PhasePartition(ts_eq).ice ≈ PhasePartition(ts_dry).ice
 
+    @test mixing_ratios(ts_eq).tot ≈ mixing_ratios(ts_dry).tot
+    @test mixing_ratios(ts_eq).liq ≈ mixing_ratios(ts_dry).liq
+    @test mixing_ratios(ts_eq).ice ≈ mixing_ratios(ts_dry).ice
+
     ts_dry = PhaseDry.(param_set, e_int, ρ)
     ts_eq = PhaseEquil.(param_set, e_int, ρ, q_tot .* 0)