diff --git a/config/model_configs/prognostic_edmfx_trmm_column_0M.yml b/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
index 5f8db48bae..2634540e9f 100644
--- a/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
+++ b/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
@@ -30,7 +30,7 @@ z_elem: 82
 z_stretch: false
 dz_bottom: 30
 dt: 150secs
-t_end: 6hours
+t_end: 12hours
 toml: [toml/prognostic_edmfx_implicit_scm_calibrated_5_cases_shallow_deep_v1.toml]
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl
index faec9fee79..f310566528 100644
--- a/reproducibility_tests/ref_counter.jl
+++ b/reproducibility_tests/ref_counter.jl
@@ -1,4 +1,4 @@
-275
+276
 
 # **README**
 #
@@ -20,6 +20,10 @@
 
 
 #=
+276
+- Update prognostic EDMF boundary conditions: apply equal surface fluxes to the
+  updraft and grid mean, and enable entrainment of buoyant air in the first cell.
+
 275
 - Change order of GPU calculations for better performance, but it
   results in slightly different floating point rounding.  Artifacts
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 0f9d00fb85..6c40f386c2 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -556,9 +556,7 @@ NVTX.@annotate function set_explicit_precomputed_quantities_part1!(Y, p, t)
             ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts)))
     end
 
-    if turbconv_model isa PrognosticEDMFX
-        set_prognostic_edmf_precomputed_quantities_bottom_bc!(Y, p, t)
-    elseif turbconv_model isa DiagnosticEDMFX
+    if turbconv_model isa DiagnosticEDMFX
         set_diagnostic_edmf_precomputed_quantities_bottom_bc!(Y, p, t)
     elseif !(isnothing(turbconv_model))
         # Do nothing for other turbconv models for now
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index 9322b93a86..cc9c1e7d03 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -115,150 +115,6 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_draft!(
     return nothing
 end
 
-"""
-    set_prognostic_edmf_precomputed_quantities_bottom_bc!(Y, p, ᶠuₕ³, t)
-
-Updates velocity and thermodynamics quantities at the surface in each SGS draft.
-"""
-NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_bottom_bc!(
-    Y,
-    p,
-    t,
-)
-    (; moisture_model, turbconv_model, microphysics_model) = p.atmos
-
-    FT = Spaces.undertype(axes(Y.c))
-    n = n_mass_flux_subdomains(turbconv_model)
-    thermo_params = CAP.thermodynamics_params(p.params)
-    turbconv_params = CAP.turbconv_params(p.params)
-
-    (; ᶜΦ,) = p.core
-    (; ᶜp, ᶜK, ᶜtsʲs, ᶜρʲs, ᶜts) = p.precomputed
-    (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
-
-    for j in 1:n
-        ᶜtsʲ = ᶜtsʲs.:($j)
-        ᶜmseʲ = Y.c.sgsʲs.:($j).mse
-        ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
-        if moisture_model isa NonEquilMoistModel && (
-            microphysics_model isa Microphysics1Moment ||
-            microphysics_model isa Microphysics2Moment
-        )
-            ᶜq_liqʲ = Y.c.sgsʲs.:($j).q_liq
-            ᶜq_iceʲ = Y.c.sgsʲs.:($j).q_ice
-            ᶜq_raiʲ = Y.c.sgsʲs.:($j).q_rai
-            ᶜq_snoʲ = Y.c.sgsʲs.:($j).q_sno
-        end
-
-        # We need field_values everywhere because we are mixing
-        # information from surface and first interior inside the
-        # sgs_scalar_first_interior_bc call.
-        ᶜz_int_val =
-            Fields.field_values(Fields.level(Fields.coordinate_field(Y.c).z, 1))
-        z_sfc_val = Fields.field_values(
-            Fields.level(Fields.coordinate_field(Y.f).z, Fields.half),
-        )
-        ᶜρ_int_val = Fields.field_values(Fields.level(Y.c.ρ, 1))
-        ᶜp_int_val = Fields.field_values(Fields.level(ᶜp, 1))
-        (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length) =
-            p.precomputed.sfc_conditions
-
-        buoyancy_flux_val = Fields.field_values(buoyancy_flux)
-        ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
-        ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
-
-        ustar_val = Fields.field_values(ustar)
-        obukhov_length_val = Fields.field_values(obukhov_length)
-        sfc_local_geometry_val = Fields.field_values(
-            Fields.local_geometry_field(Fields.level(Y.f, Fields.half)),
-        )
-
-        # Based on boundary conditions for updrafts we overwrite
-        # the first interior point for EDMFX ᶜmseʲ...
-        ᶜaʲ_int_val = p.scratch.temp_data_level
-        # TODO: replace this with the actual surface area fraction when
-        # using prognostic surface area
-        @. ᶜaʲ_int_val = FT(turbconv_params.surface_area)
-        ᶜh_tot = @. lazy(
-            TD.total_specific_enthalpy(
-                thermo_params,
-                ᶜts,
-                specific(Y.c.ρe_tot, Y.c.ρ),
-            ),
-        )
-        ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
-        ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1))
-        ᶜmseʲ_int_val = Fields.field_values(Fields.level(ᶜmseʲ, 1))
-        @. ᶜmseʲ_int_val = sgs_scalar_first_interior_bc(
-            ᶜz_int_val - z_sfc_val,
-            ᶜρ_int_val,
-            ᶜaʲ_int_val,
-            ᶜh_tot_int_val - ᶜK_int_val,
-            buoyancy_flux_val,
-            ρ_flux_h_tot_val,
-            ustar_val,
-            obukhov_length_val,
-            sfc_local_geometry_val,
-        )
-
-        # ... and the first interior point for EDMFX ᶜq_totʲ.
-
-        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
-        ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
-        ᶜq_totʲ_int_val = Fields.field_values(Fields.level(ᶜq_totʲ, 1))
-        @. ᶜq_totʲ_int_val = sgs_scalar_first_interior_bc(
-            ᶜz_int_val - z_sfc_val,
-            ᶜρ_int_val,
-            ᶜaʲ_int_val,
-            ᶜq_tot_int_val,
-            buoyancy_flux_val,
-            ρ_flux_q_tot_val,
-            ustar_val,
-            obukhov_length_val,
-            sfc_local_geometry_val,
-        )
-
-        # Then overwrite the prognostic variables at first inetrior point.
-        ᶜΦ_int_val = Fields.field_values(Fields.level(ᶜΦ, 1))
-        ᶜtsʲ_int_val = Fields.field_values(Fields.level(ᶜtsʲ, 1))
-        if moisture_model isa NonEquilMoistModel && (
-            microphysics_model isa Microphysics1Moment ||
-            microphysics_model isa Microphysics2Moment
-        )
-            ᶜq_liqʲ_int_val = Fields.field_values(Fields.level(ᶜq_liqʲ, 1))
-            ᶜq_iceʲ_int_val = Fields.field_values(Fields.level(ᶜq_iceʲ, 1))
-            ᶜq_raiʲ_int_val = Fields.field_values(Fields.level(ᶜq_raiʲ, 1))
-            ᶜq_snoʲ_int_val = Fields.field_values(Fields.level(ᶜq_snoʲ, 1))
-            @. ᶜtsʲ_int_val = TD.PhaseNonEquil_phq(
-                thermo_params,
-                ᶜp_int_val,
-                ᶜmseʲ_int_val - ᶜΦ_int_val,
-                TD.PhasePartition(
-                    ᶜq_totʲ_int_val,
-                    ᶜq_liqʲ_int_val + ᶜq_raiʲ_int_val,
-                    ᶜq_iceʲ_int_val + ᶜq_snoʲ_int_val,
-                ),
-            )
-        else
-            @. ᶜtsʲ_int_val = TD.PhaseEquil_phq(
-                thermo_params,
-                ᶜp_int_val,
-                ᶜmseʲ_int_val - ᶜΦ_int_val,
-                ᶜq_totʲ_int_val,
-            )
-        end
-        sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1))
-        sgsʲs_ρa_int_val =
-            Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1))
-
-        @. sgsʲs_ρ_int_val = TD.air_density(thermo_params, ᶜtsʲ_int_val)
-        @. sgsʲs_ρa_int_val =
-            $(FT(turbconv_params.surface_area)) *
-            TD.air_density(thermo_params, ᶜtsʲ_int_val)
-    end
-    return nothing
-end
-
 """
     set_prognostic_edmf_precomputed_quantities_implicit_closures!(Y, p, t)
 
@@ -332,6 +188,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
         ᶠnh_pressure₃_buoyʲs,
     ) = p.precomputed
     (; ustar, obukhov_length) = p.precomputed.sfc_conditions
+    ᶜaʲ_int_val = p.scratch.temp_data_level
 
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
@@ -411,6 +268,27 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos
             dt,
         )
 
+        # If the surface buoyancy flux is positive, increase entrainment in the first cell
+        # so that the updraft area grows to at least `surface_area` within one timestep.
+        # Otherwise (stable surface), leave entrainment unchanged.
+        buoyancy_flux_val = Fields.field_values(p.precomputed.sfc_conditions.buoyancy_flux)
+        sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1))
+        sgsʲs_ρa_int_val =
+            Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1))
+        @. ᶜaʲ_int_val = draft_area(sgsʲs_ρa_int_val, sgsʲs_ρ_int_val)
+        entr_int_val = Fields.field_values(Fields.level(ᶜentrʲs.:($j), 1))
+        detr_int_val = Fields.field_values(Fields.level(ᶜdetrʲs.:($j), 1))
+        @. entr_int_val = limit_entrainment(
+            ifelse(
+                buoyancy_flux_val < 0 || ᶜaʲ_int_val >= $(FT(turbconv_params.surface_area)),
+                entr_int_val,
+                detr_int_val +
+                ($(FT(turbconv_params.surface_area)) / ᶜaʲ_int_val - 1) / FT(dt),
+            ),
+            ᶜaʲ_int_val,
+            dt,
+        )
+
         # The buoyancy term in the nonhydrostatic pressure closure is always applied
         # for prognostic edmf. The tendency is combined with the buoyancy term in the
         # updraft momentum equation in `edmfx_sgs_vertical_advection_tendency!`. This
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index be3b214bb1..b110c64404 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -580,6 +580,139 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMF
     return nothing
 end
 
+"""
+    edmfx_first_interior_entr_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX)
+
+Apply first-interior–level entrainment tendencies for each EDMF updraft.
+
+This routine (1) seeds a small positive updraft area fraction when surface
+buoyancy flux is positive—allowing the plume to grow from zero—and  
+(2) adds entrainment tendencies for moist static energy (`mse`) and total
+humidity (`q_tot`) in the first model cell.  
+The entrained tracer value is taken from `sgs_scalar_first_interior_bc`.
+"""
+edmfx_first_interior_entr_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
+function edmfx_first_interior_entr_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    turbconv_model::PrognosticEDMFX,
+)
+
+    (; params, dt) = p
+    (; ᶜK, ᶜρʲs, ᶜentrʲs, ᶜts) = p.precomputed
+
+    FT = eltype(params)
+    n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+    thermo_params = CAP.thermodynamics_params(params)
+    turbconv_params = CAP.turbconv_params(params)
+    ᶜaʲ_int_val = p.scratch.temp_data_level
+    ᶜmse_buoyant_air_int_val = p.scratch.temp_data_level_2
+    ᶜq_tot_buoyant_air_int_val = p.scratch.temp_data_level_3
+
+    (;
+        ustar,
+        obukhov_length,
+        buoyancy_flux,
+        ρ_flux_h_tot,
+        ρ_flux_q_tot,
+        ustar,
+        obukhov_length,
+    ) =
+        p.precomputed.sfc_conditions
+
+    ᶜz_int_val = Fields.field_values(Fields.level(Fields.coordinate_field(Y.c).z, 1))
+    z_sfc_val =
+        Fields.field_values(Fields.level(Fields.coordinate_field(Y.f).z, Fields.half))
+    ᶜρ_int_val = Fields.field_values(Fields.level(Y.c.ρ, 1))
+
+    buoyancy_flux_val = Fields.field_values(buoyancy_flux)
+    ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
+    ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
+
+    ustar_val = Fields.field_values(ustar)
+    obukhov_length_val = Fields.field_values(obukhov_length)
+    sfc_local_geometry_val = Fields.field_values(
+        Fields.local_geometry_field(Fields.level(Y.f, Fields.half)),
+    )
+
+    ᶜh_tot = @. lazy(
+        TD.total_specific_enthalpy(
+            thermo_params,
+            ᶜts,
+            specific(Y.c.ρe_tot, Y.c.ρ),
+        ),
+    )
+    ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
+    ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1))
+    ᶜmse⁰ = ᶜspecific_env_mse(Y, p)
+    env_mse_int_val = Fields.field_values(Fields.level(ᶜmse⁰, 1))
+
+    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+    ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
+    ᶜq_tot⁰ = ᶜspecific_env_value(@name(q_tot), Y, p)
+    env_q_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot⁰, 1))
+
+    for j in 1:n
+
+        # Seed a small positive updraft area fraction when surface buoyancy flux is positive.
+        # This perturbation prevents the plume area from staying identically zero,
+        # allowing entrainment to grow it to the prescribed surface area.
+        sgsʲs_ρ_int_val = Fields.field_values(Fields.level(ᶜρʲs.:($j), 1))
+        sgsʲs_ρa_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1))
+        sgsʲs_ρaₜ_int_val = Fields.field_values(Fields.level(Yₜ.c.sgsʲs.:($j).ρa, 1))
+        @. sgsʲs_ρaₜ_int_val += ifelse(buoyancy_flux_val < 0,
+            0,
+            max(0, (sgsʲs_ρ_int_val * $(eps(FT)) - sgsʲs_ρa_int_val) / FT(dt)),
+        )
+
+        # Apply entrainment tendencies in the first model cell for moist static energy (mse) 
+        # and total humidity (q_tot). The entrained fluid is assumed to have a scalar value 
+        # given by `sgs_scalar_first_interior_bc` (mean + SGS perturbation). Since 
+        # `edmfx_entr_detr_tendency!` computes entrainment based on the environment–updraft 
+        # contrast, we supply the high-value (entrained) tracer minus the environment value 
+        # here to form the correct tendency.
+        entr_int_val = Fields.field_values(Fields.level(ᶜentrʲs.:($j), 1))
+        @. ᶜaʲ_int_val = max(
+            FT(turbconv_params.surface_area),
+            draft_area(sgsʲs_ρa_int_val, sgsʲs_ρ_int_val),
+        )
+
+        sgsʲs_mseₜ_int_val =
+            Fields.field_values(Fields.level(Yₜ.c.sgsʲs.:($j).mse, 1))
+        @. ᶜmse_buoyant_air_int_val = sgs_scalar_first_interior_bc(
+            ᶜz_int_val - z_sfc_val,
+            ᶜρ_int_val,
+            ᶜaʲ_int_val,
+            ᶜh_tot_int_val - ᶜK_int_val,
+            buoyancy_flux_val,
+            ρ_flux_h_tot_val,
+            ustar_val,
+            obukhov_length_val,
+            sfc_local_geometry_val,
+        )
+        @. sgsʲs_mseₜ_int_val += entr_int_val * (ᶜmse_buoyant_air_int_val - env_mse_int_val)
+
+        sgsʲs_q_totₜ_int_val =
+            Fields.field_values(Fields.level(Yₜ.c.sgsʲs.:($j).q_tot, 1))
+        @. ᶜq_tot_buoyant_air_int_val = sgs_scalar_first_interior_bc(
+            ᶜz_int_val - z_sfc_val,
+            ᶜρ_int_val,
+            ᶜaʲ_int_val,
+            ᶜq_tot_int_val,
+            buoyancy_flux_val,
+            ρ_flux_q_tot_val,
+            ustar_val,
+            obukhov_length_val,
+            sfc_local_geometry_val,
+        )
+        @. sgsʲs_q_totₜ_int_val +=
+            entr_int_val * (ᶜq_tot_buoyant_air_int_val - env_q_tot_int_val)
+
+    end
+end
+
 # limit entrainment and detrainment rates for prognostic EDMFX
 # limit rates approximately below the inverse timescale 1/dt
 limit_entrainment(entr::FT, a, dt) where {FT} = max(
diff --git a/src/prognostic_equations/mass_flux_closures.jl b/src/prognostic_equations/mass_flux_closures.jl
index 9a5a4250fd..8026644b5c 100644
--- a/src/prognostic_equations/mass_flux_closures.jl
+++ b/src/prognostic_equations/mass_flux_closures.jl
@@ -206,6 +206,7 @@ end
 Apply EDMF filters:
  - Relax u_3 to zero when it is negative
  - Relax ρa to zero when it is negative
+ - Relax ρa to ρ when a is larger than one
 
 This function modifies the tendency `Yₜ` in place based on the current state `Y`,
 parameters `p`, time `t`, and the turbulence convection model `turbconv_model`.
@@ -225,6 +226,8 @@ function edmfx_filter_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX)
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
                 C3(min(Y.f.sgsʲs.:($$j).u₃.components.data.:1, 0)) / FT(dt)
             @. Yₜ.c.sgsʲs.:($$j).ρa -= min(Y.c.sgsʲs.:($$j).ρa, 0) / FT(dt)
+            @. Yₜ.c.sgsʲs.:($$j).ρa -=
+                max(Y.c.sgsʲs.:($$j).ρa - p.precomputed.ᶜρʲs.:($$j), 0) / FT(dt)
         end
     end
 end
diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
index 9160fe34d6..3428df15eb 100644
--- a/src/prognostic_equations/remaining_tendency.jl
+++ b/src/prognostic_equations/remaining_tendency.jl
@@ -284,6 +284,7 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     if p.atmos.sgs_entr_detr_mode == Explicit()
         edmfx_entr_detr_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
+    edmfx_first_interior_entr_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     if p.atmos.sgs_mf_mode == Explicit()
         edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
diff --git a/src/prognostic_equations/surface_flux.jl b/src/prognostic_equations/surface_flux.jl
index e76e083040..e42ccccee3 100644
--- a/src/prognostic_equations/surface_flux.jl
+++ b/src/prognostic_equations/surface_flux.jl
@@ -105,6 +105,7 @@ function surface_flux_tendency!(Yₜ, Y, p, t)
 
     FT = eltype(Y)
     (; params) = p
+    (; turbconv_model) = p.atmos
     (; sfc_conditions, ᶜts) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
 
@@ -123,6 +124,12 @@ function surface_flux_tendency!(Yₜ, Y, p, t)
     )
     btt = boundary_tendency_scalar(ᶜh_tot, sfc_conditions.ρ_flux_h_tot)
     @. Yₜ.c.ρe_tot -= btt
+
+    if turbconv_model isa PrognosticEDMFX
+        # assuming one updraft
+        @. Yₜ.c.sgsʲs.:(1).mse -= specific(btt, p.precomputed.ᶜρʲs.:(1))
+    end
+
     ρ_flux_χ = p.scratch.sfc_temp_C3
     foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name
         ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ))
@@ -136,5 +143,11 @@ function surface_flux_tendency!(Yₜ, Y, p, t)
         if ρχ_name == @name(ρq_tot)
             @. Yₜ.c.ρ -= btt
         end
+
+        if turbconv_model isa PrognosticEDMFX
+            # assuming one updraft
+            ᶜχʲₜ = MatrixFields.get_field(Yₜ.c, get_χʲ_name_from_ρχ_name(ρχ_name))
+            @. ᶜχʲₜ -= specific(btt, p.precomputed.ᶜρʲs.:(1))
+        end
     end
 end