diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index f142189db1..57e48cd631 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -175,6 +175,9 @@ divergence_damping_factor:
 rayleigh_sponge:
   help: "Rayleigh sponge [`true`, `false` (default)]"
   value: false
+disable_surface_flux_tendency:
+  help: "(Bool) Whether to disable surface flux tendencies of momentum, energy, and tracers [`true`, `false` (default)]. When this flag is true, the surface flux tendency is not applied, no matter how surface conditions are computed."
+  value: false
 surface_setup:
   help: "Surface flux scheme [`DefaultExchangeCoefficients` (default), `DefaultMoninObukhov`]"
   value: "DefaultExchangeCoefficients"
diff --git a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml
index 78f4d558c4..f9b274bcc7 100644
--- a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml
+++ b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean.yml
@@ -11,4 +11,4 @@ rad: "clearsky"
 insolation: "timevarying"
 dt_rad: "1hours"
 prognostic_surface: "PrognosticSurfaceTemperature"
-#check_conservation: true
+check_conservation: true
diff --git a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml
index b70a5dafc1..10be964983 100644
--- a/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml
+++ b/config/model_configs/aquaplanet_equil_clearsky_tvinsol_0M_slabocean_ft64.yml
@@ -11,4 +11,4 @@ rad: "clearsky"
 insolation: "timevarying"
 dt_rad: "1hours"
 prognostic_surface: "PrognosticSurfaceTemperature"
-#check_conservation: true
+check_conservation: true
diff --git a/config/model_configs/baroclinic_wave.yml b/config/model_configs/baroclinic_wave.yml
index 210c39b46c..0e36935273 100644
--- a/config/model_configs/baroclinic_wave.yml
+++ b/config/model_configs/baroclinic_wave.yml
@@ -3,6 +3,7 @@ initial_condition: "DryBaroclinicWave"
 dt: "400secs"
 t_end: "10days"
 deep_atmosphere: false
+disable_surface_flux_tendency: true
 diagnostics:
   - short_name: [pfull, ua, wa, va, rv, ta, ke]
     period: 1days
diff --git a/config/model_configs/baroclinic_wave_conservation.yml b/config/model_configs/baroclinic_wave_conservation.yml
index 29c6794d7a..974abb85fc 100644
--- a/config/model_configs/baroclinic_wave_conservation.yml
+++ b/config/model_configs/baroclinic_wave_conservation.yml
@@ -2,6 +2,7 @@ initial_condition: "DryBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 dt_save_state_to_disk: "5days"
+disable_surface_flux_tendency: true
 check_conservation: true
 diagnostics:
   - short_name: [massa, energya]
diff --git a/config/model_configs/baroclinic_wave_equil.yml b/config/model_configs/baroclinic_wave_equil.yml
index 5e5079aa77..dbd4527d23 100644
--- a/config/model_configs/baroclinic_wave_equil.yml
+++ b/config/model_configs/baroclinic_wave_equil.yml
@@ -5,6 +5,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "450secs"
 t_end: "10days"
 moist: "equil"
+disable_surface_flux_tendency: true
 diagnostics:
   - short_name: [pfull, wa, va, rv, hus, ke]
     period: 1days
diff --git a/config/model_configs/baroclinic_wave_equil_conservation.yml b/config/model_configs/baroclinic_wave_equil_conservation.yml
index 4636d419a8..261a20ed9e 100644
--- a/config/model_configs/baroclinic_wave_equil_conservation.yml
+++ b/config/model_configs/baroclinic_wave_equil_conservation.yml
@@ -2,6 +2,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 moist: "equil"
+disable_surface_flux_tendency: true
 dt_save_state_to_disk: "5days"
 check_conservation: true
 diagnostics:
diff --git a/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml b/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml
index 2d2440304b..310d7eedb7 100644
--- a/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml
+++ b/config/model_configs/baroclinic_wave_equil_conservation_ft64.yml
@@ -3,6 +3,7 @@ initial_condition: "MoistBaroclinicWave"
 dt: "400secs"
 t_end: "5days"
 moist: "equil"
+disable_surface_flux_tendency: true
 dt_save_state_to_disk: "5days"
 check_conservation: true
 diagnostics:
diff --git a/config/model_configs/baroclinic_wave_hughes2023.yml b/config/model_configs/baroclinic_wave_hughes2023.yml
index 64d407af43..9c2c710a8a 100644
--- a/config/model_configs/baroclinic_wave_hughes2023.yml
+++ b/config/model_configs/baroclinic_wave_hughes2023.yml
@@ -3,6 +3,7 @@ initial_condition: "MoistBaroclinicWave"
 topography: "Hughes2023"
 topo_smoothing: false
 moist: equil
+disable_surface_flux_tendency: true
 dt: "200secs"
 t_end: "10days"
 diagnostics:
diff --git a/config/model_configs/baroclinic_wave_topography_dcmip_rs.yml b/config/model_configs/baroclinic_wave_topography_dcmip_rs.yml
index 4563d5af27..ff3eec2ad9 100644
--- a/config/model_configs/baroclinic_wave_topography_dcmip_rs.yml
+++ b/config/model_configs/baroclinic_wave_topography_dcmip_rs.yml
@@ -6,6 +6,7 @@ initial_condition: "DryBaroclinicWave"
 t_end: "6days"
 dt: "200secs"
 hyperdiff: "ClimaHyperdiffusion"
+disable_surface_flux_tendency: true
 netcdf_output_at_levels: false
 diagnostics:
   - short_name: [pfull, wa, va, rv]
diff --git a/config/model_configs/box_cosine_hills_float64_test.yml b/config/model_configs/box_cosine_hills_float64_test.yml
index 1b4948a974..739a6a05bb 100644
--- a/config/model_configs/box_cosine_hills_float64_test.yml
+++ b/config/model_configs/box_cosine_hills_float64_test.yml
@@ -12,6 +12,7 @@ dz_bottom: 10
 dt: "0.5secs" # CFL is slightly lower than for 2D case
 t_end: "30mins"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/box_density_current_test.yml b/config/model_configs/box_density_current_test.yml
index 11805262c2..42fdf89e49 100644
--- a/config/model_configs/box_density_current_test.yml
+++ b/config/model_configs/box_density_current_test.yml
@@ -14,6 +14,7 @@ z_stretch: false
 dt: "0.25secs"
 t_end: "20secs"
 dt_save_state_to_disk: "5secs"
+disable_surface_flux_tendency: true
 netcdf_interpolation_num_points: [40, 40, 80]
 diagnostics:
   - short_name: [wa, ua, va, ta, thetaa, ha]
diff --git a/config/model_configs/box_hydrostatic_balance.yml b/config/model_configs/box_hydrostatic_balance.yml
index a4d2bebfff..63ecc51d06 100644
--- a/config/model_configs/box_hydrostatic_balance.yml
+++ b/config/model_configs/box_hydrostatic_balance.yml
@@ -3,6 +3,7 @@ config: "box"
 hyperdiff: false
 dt: "20secs"
 perturb_initstate: false
+disable_surface_flux_tendency: true
 diagnostics:
   - short_name: [wa, ua]
     period: 1days
diff --git a/config/model_configs/extruded_plane_cosine_hills_float64_test.yml b/config/model_configs/extruded_plane_cosine_hills_float64_test.yml
index 8542558ddb..100e15bd9a 100644
--- a/config/model_configs/extruded_plane_cosine_hills_float64_test.yml
+++ b/config/model_configs/extruded_plane_cosine_hills_float64_test.yml
@@ -12,6 +12,7 @@ dz_bottom: 10
 dt: "0.5secs" # CFL is slightly lower than for 2D case
 t_end: "9hours"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/hydrostatic_balance_ft64.yml b/config/model_configs/hydrostatic_balance_ft64.yml
index f5c48d0c37..cd9f33efef 100644
--- a/config/model_configs/hydrostatic_balance_ft64.yml
+++ b/config/model_configs/hydrostatic_balance_ft64.yml
@@ -4,4 +4,5 @@ t_end: "8days"
 discrete_hydrostatic_balance: true
 hyperdiff: false
 perturb_initstate: false
+disable_surface_flux_tendency: true
 FLOAT_TYPE: "Float64"
diff --git a/config/model_configs/plane_agnesi_mountain_float64_test.yml b/config/model_configs/plane_agnesi_mountain_float64_test.yml
index 0ab10e94e7..e84a19249b 100644
--- a/config/model_configs/plane_agnesi_mountain_float64_test.yml
+++ b/config/model_configs/plane_agnesi_mountain_float64_test.yml
@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/plane_cosine_hills_float64_test.yml b/config/model_configs/plane_cosine_hills_float64_test.yml
index d5b8ef67af..b50c71c709 100644
--- a/config/model_configs/plane_cosine_hills_float64_test.yml
+++ b/config/model_configs/plane_cosine_hills_float64_test.yml
@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/plane_density_current_test.yml b/config/model_configs/plane_density_current_test.yml
index a982d9699c..24c2cb4bcb 100644
--- a/config/model_configs/plane_density_current_test.yml
+++ b/config/model_configs/plane_density_current_test.yml
@@ -9,6 +9,7 @@ z_stretch: false
 x_elem: 80
 config: "plane"
 z_max: 6400.0
+disable_surface_flux_tendency: true
 output_default_diagnostics: false
 diagnostics:
   - short_name: thetaa
diff --git a/config/model_configs/plane_no_topography_float64_test.yml b/config/model_configs/plane_no_topography_float64_test.yml
index 7cd38353cd..5458922c05 100644
--- a/config/model_configs/plane_no_topography_float64_test.yml
+++ b/config/model_configs/plane_no_topography_float64_test.yml
@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/plane_schar_mountain_float32_test.yml b/config/model_configs/plane_schar_mountain_float32_test.yml
index 99799a5488..3532336fe4 100644
--- a/config/model_configs/plane_schar_mountain_float32_test.yml
+++ b/config/model_configs/plane_schar_mountain_float32_test.yml
@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/plane_schar_mountain_float64_test.yml b/config/model_configs/plane_schar_mountain_float64_test.yml
index 9f912768c5..e62fb575ea 100644
--- a/config/model_configs/plane_schar_mountain_float64_test.yml
+++ b/config/model_configs/plane_schar_mountain_float64_test.yml
@@ -10,6 +10,7 @@ dz_bottom: 10
 dt: "0.7secs"
 t_end: "2days"
 rayleigh_sponge: true
+disable_surface_flux_tendency: true
 toml: [toml/steady_state_test.toml]
 check_steady_state: true
 output_default_diagnostics: false
diff --git a/config/model_configs/single_column_hydrostatic_balance_ft64.yml b/config/model_configs/single_column_hydrostatic_balance_ft64.yml
index 549ecf4a69..cea418d9d1 100644
--- a/config/model_configs/single_column_hydrostatic_balance_ft64.yml
+++ b/config/model_configs/single_column_hydrostatic_balance_ft64.yml
@@ -5,6 +5,7 @@ dt_save_state_to_disk: "5days"
 initial_condition: "IsothermalProfile"
 config: "column"
 hyperdiff: false
+disable_surface_flux_tendency: true
 dt: "3hours"
 FLOAT_TYPE: "Float64"
 reproducibility_test: true
diff --git a/config/model_configs/single_column_nonorographic_gravity_wave.yml b/config/model_configs/single_column_nonorographic_gravity_wave.yml
index d36e27440a..de44e5d498 100644
--- a/config/model_configs/single_column_nonorographic_gravity_wave.yml
+++ b/config/model_configs/single_column_nonorographic_gravity_wave.yml
@@ -4,4 +4,5 @@ t_end: "1500secs"
 config: "column"
 hyperdiff: false
 dt: "400secs"
+disable_surface_flux_tendency: true
 non_orographic_gravity_wave: true
diff --git a/config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml b/config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml
index 80e7b6f14c..0e262a0353 100644
--- a/config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml
+++ b/config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml
@@ -12,6 +12,7 @@ idealized_clouds: true
 config: "column"
 insolation: "timevarying"
 dt_save_to_sol: "30hours"
+disable_surface_flux_tendency: true
 rad: "allskywithclear"
 diagnostics:
   - short_name: rhoa
diff --git a/config/model_configs/single_column_radiative_equilibrium_clearsky.yml b/config/model_configs/single_column_radiative_equilibrium_clearsky.yml
index d8476317b5..a2b97e6f2d 100644
--- a/config/model_configs/single_column_radiative_equilibrium_clearsky.yml
+++ b/config/model_configs/single_column_radiative_equilibrium_clearsky.yml
@@ -10,4 +10,5 @@ idealized_h2o: true
 t_end: "654days"
 config: "column"
 dt_save_to_sol: "30hours"
+disable_surface_flux_tendency: true
 rad: "clearsky"
diff --git a/config/model_configs/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.yml b/config/model_configs/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.yml
index e9bf2619b4..bbb05611b5 100644
--- a/config/model_configs/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.yml
+++ b/config/model_configs/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.yml
@@ -11,5 +11,6 @@ dt: "3hours"
 dt_rad: "3hours"
 dt_save_to_sol: "30hours"
 dt_save_state_to_disk: "100days"
+disable_surface_flux_tendency: true
 prognostic_surface: true
 toml: [toml/single_column_radiative_equilibrium_clearsky_prognostic_surface_temp.toml]
diff --git a/config/model_configs/single_column_radiative_equilibrium_gray.yml b/config/model_configs/single_column_radiative_equilibrium_gray.yml
index e3ddb53a04..256f039fda 100644
--- a/config/model_configs/single_column_radiative_equilibrium_gray.yml
+++ b/config/model_configs/single_column_radiative_equilibrium_gray.yml
@@ -10,6 +10,7 @@ t_end: "654days"
 dz_bottom: 30.0
 config: "column"
 dt_save_to_sol: "30hours"
+disable_surface_flux_tendency: true
 rad: "gray"
 # [2, 2, 80] instead of [1, 1, 80] because Julia ranges are inclusive of the
 # extrema. Given that our columns are 3D, we cannot map the horizontal dimension
diff --git a/config/model_configs/sphere_nonorographic_gravity_wave.yml b/config/model_configs/sphere_nonorographic_gravity_wave.yml
index 89b93e0408..bdb8c04298 100644
--- a/config/model_configs/sphere_nonorographic_gravity_wave.yml
+++ b/config/model_configs/sphere_nonorographic_gravity_wave.yml
@@ -1,3 +1,4 @@
 t_end: "1500secs"
 dt: "400secs"
+disable_surface_flux_tendency: true
 non_orographic_gravity_wave: true
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
index 46d739cc26..05d727ca12 100644
--- a/src/ClimaAtmos.jl
+++ b/src/ClimaAtmos.jl
@@ -116,6 +116,7 @@ include(
 include(
     joinpath("prognostic_equations", "vertical_diffusion_boundary_layer.jl"),
 )
+include(joinpath("prognostic_equations", "surface_flux.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "rayleigh_sponge.jl"))
 include(joinpath("parameterized_tendencies", "sponge", "viscous_sponge.jl"))
 include(
diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
index c36fd6b070..b1dceb623f 100644
--- a/src/parameters/create_parameters.jl
+++ b/src/parameters/create_parameters.jl
@@ -36,7 +36,7 @@ function ClimaAtmosParameters(toml_dict::TD) where {TD <: CP.AbstractTOMLDict}
     IP = typeof(insolation_params)
 
     surface_fluxes_params =
-        SF.Parameters.SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
+        SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
     SFP = typeof(surface_fluxes_params)
 
     surface_temp_params = SurfaceTemperatureParameters(toml_dict)
diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
index 850450854e..2f479b2cbf 100644
--- a/src/prognostic_equations/edmfx_sgs_flux.jl
+++ b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -178,7 +178,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; sfc_conditions) = p.precomputed
     (; ᶜρa⁰, ᶜu⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
@@ -192,7 +191,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
         # energy
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜmse⁰ + ᶜK⁰)))
         if use_prognostic_tke(turbconv_model)
@@ -216,7 +215,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
             ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
             ᶜdivᵥ_ρq_tot = Operators.DivergenceF2C(
                 top = Operators.SetValue(C3(FT(0))),
-                bottom = Operators.SetValue(sfc_conditions.ρ_flux_q_tot),
+                bottom = Operators.SetValue(C3(FT(0))),
             )
             @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜq_tot⁰)))
             @. Yₜ.c.ρq_tot -= ᶜρχₜ_diffusion
@@ -228,12 +227,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
         bc_strain_rate = compute_strain_rate_face(ᶜu⁰)
         @. ᶠstrain_rate = bc_strain_rate
         @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(-(2 * ᶠρaK_u * ᶠstrain_rate)) / Y.c.ρ)
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     # TODO: Add tracer flux
@@ -253,7 +246,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
     (; dt, params) = p
     turbconv_params = CAP.turbconv_params(params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    (; sfc_conditions) = p.precomputed
     (; ᶜu, ᶜh_tot, ᶜspecific, ᶜtke⁰, ᶜmixing_length) = p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
     ᶠgradᵥ = Operators.GradientC2F()
@@ -267,7 +259,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
         # energy
         ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜh_tot)))
 
@@ -288,7 +280,7 @@ function edmfx_sgs_diffusive_flux_tendency!(
             ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
             ᶜdivᵥ_ρq_tot = Operators.DivergenceF2C(
                 top = Operators.SetValue(C3(FT(0))),
-                bottom = Operators.SetValue(sfc_conditions.ρ_flux_q_tot),
+                bottom = Operators.SetValue(C3(FT(0))),
             )
             @. ᶜρχₜ_diffusion =
                 ᶜdivᵥ_ρq_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜspecific.q_tot)))
@@ -301,12 +293,6 @@ function edmfx_sgs_diffusive_flux_tendency!(
         bc_strain_rate = compute_strain_rate_face(ᶜu)
         @. ᶠstrain_rate = bc_strain_rate
         @. Yₜ.c.uₕ -= C12(ᶜdivᵥ(-(2 * ᶠρaK_u * ᶠstrain_rate)) / Y.c.ρ)
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     # TODO: Add tracer flux
diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl
index 5ba0b2d71e..039af2b64e 100644
--- a/src/prognostic_equations/remaining_tendency.jl
+++ b/src/prognostic_equations/remaining_tendency.jl
@@ -109,6 +109,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
         edmfx_sgs_diffusive_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
 
+    surface_flux_tendency!(Yₜ, Y, p, t)
+
     radiation_tendency!(Yₜ, Y, p, t, p.atmos.radiation_mode)
     edmfx_entr_detr_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     if p.atmos.sgs_mf_mode == Explicit()
diff --git a/src/prognostic_equations/surface_flux.jl b/src/prognostic_equations/surface_flux.jl
new file mode 100644
index 0000000000..f65f4b1c7a
--- /dev/null
+++ b/src/prognostic_equations/surface_flux.jl
@@ -0,0 +1,49 @@
+#####
+##### Surface flux tendencies
+#####
+
+import ClimaCore.Geometry: ⊗
+import ClimaCore.Operators as Operators
+
+function surface_flux_tendency!(Yₜ, Y, p, t)
+
+    p.atmos.disable_surface_flux_tendency && return
+    
+    FT = eltype(Y)
+    (; ᶜh_tot, ᶜspecific, sfc_conditions) = p.precomputed
+    ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ
+
+    if diffuse_momentum(p.atmos.vert_diff) || !isnothing(p.atmos.turbconv_model)
+        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
+            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
+        )
+        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
+    end
+
+    ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
+        top = Operators.SetValue(C3(FT(0))),
+        bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+    )
+    @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(FT(0) * ᶠgradᵥ(ᶜh_tot)))
+
+    ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
+    ρ_flux_χ = p.scratch.sfc_temp_C3
+    for (ᶜρχₜ, ᶜχ, χ_name) in matching_subfields(Yₜ.c, ᶜspecific)
+        χ_name == :e_tot && continue
+        if χ_name == :q_tot
+            @. ρ_flux_χ = sfc_conditions.ρ_flux_q_tot
+        else
+            @. ρ_flux_χ = C3(FT(0))
+        end
+        ᶜdivᵥ_ρχ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0))),
+            bottom = Operators.SetValue(ρ_flux_χ),
+        )
+        @. ᶜρχₜ_diffusion = ᶜdivᵥ_ρχ(-(FT(0) * ᶠgradᵥ(ᶜχ)))
+        @. ᶜρχₜ -= ᶜρχₜ_diffusion
+        if !(χ_name in (:q_rai, :q_sno))
+            @. Yₜ.c.ρ -= ᶜρχₜ_diffusion
+        end
+    end
+end
diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
index 0ca9a49b5a..0cf321ec78 100644
--- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
+++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl
@@ -2,15 +2,7 @@
 ##### Vertical diffusion boundary layer parameterization
 #####
 
-import StaticArrays
 import ClimaCore.Geometry: ⊗
-import ClimaCore.Utilities: half
-import LinearAlgebra: norm
-import Thermodynamics as TD
-import SurfaceFluxes as SF
-import ClimaCore.Spaces as Spaces
-import ClimaCore.Geometry as Geometry
-import ClimaCore.Fields as Fields
 import ClimaCore.Operators as Operators
 
 vertical_diffusion_boundary_layer_tendency!(Yₜ, Y, p, t) =
@@ -36,34 +28,21 @@ function vertical_diffusion_boundary_layer_tendency!(
         @. Yₜ.c.uₕ -= C12(
             ᶜdivᵥ(-2 * ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_u) * ᶠstrain_rate) / Y.c.ρ,
         )
-
-        # apply boundary condition for momentum flux
-        ᶜdivᵥ_uₕ = Operators.DivergenceF2C(
-            top = Operators.SetValue(C3(FT(0)) ⊗ C12(FT(0), FT(0))),
-            bottom = Operators.SetValue(sfc_conditions.ρ_flux_uₕ),
-        )
-        @. Yₜ.c.uₕ -= ᶜdivᵥ_uₕ(-(FT(0) * ᶠgradᵥ(Y.c.uₕ))) / Y.c.ρ
     end
 
     ᶜdivᵥ_ρe_tot = Operators.DivergenceF2C(
         top = Operators.SetValue(C3(FT(0))),
-        bottom = Operators.SetValue(sfc_conditions.ρ_flux_h_tot),
+        bottom = Operators.SetValue(C3(FT(0))),
     )
     @. Yₜ.c.ρe_tot -=
         ᶜdivᵥ_ρe_tot(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜh_tot)))
 
     ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar
-    ρ_flux_χ = p.scratch.sfc_temp_C3
     for (ᶜρχₜ, ᶜχ, χ_name) in matching_subfields(Yₜ.c, ᶜspecific)
         χ_name == :e_tot && continue
-        if χ_name == :q_tot
-            @. ρ_flux_χ = sfc_conditions.ρ_flux_q_tot
-        else
-            @. ρ_flux_χ = C3(FT(0))
-        end
         ᶜdivᵥ_ρχ = Operators.DivergenceF2C(
             top = Operators.SetValue(C3(FT(0))),
-            bottom = Operators.SetValue(ρ_flux_χ),
+            bottom = Operators.SetValue(C3(FT(0))),
         )
         @. ᶜρχₜ_diffusion =
             ᶜdivᵥ_ρχ(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜχ)))
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 86b723ce56..928d18994c 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -105,6 +105,7 @@ function get_atmos(config::AtmosConfig, params)
         rayleigh_sponge = get_rayleigh_sponge_model(parsed_args, params, FT),
         sfc_temperature = get_sfc_temperature_form(parsed_args),
         insolation = get_insolation_form(parsed_args),
+        disable_surface_flux_tendency = parsed_args["disable_surface_flux_tendency"],
         surface_model = get_surface_model(parsed_args),
         surface_albedo = get_surface_albedo_model(parsed_args, params, FT),
         numerics = get_numerics(parsed_args),
diff --git a/src/solver/types.jl b/src/solver/types.jl
index b8c9fac13d..243e89adfd 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -219,6 +219,8 @@ end
 diffuse_momentum(::DecayWithHeightDiffusion{DM}) where {DM} = DM
 diffuse_momentum(::Nothing) = false
 
+struct SurfaceFlux end
+
 abstract type AbstractSponge end
 Base.Broadcast.broadcastable(x::AbstractSponge) = tuple(x)
 Base.@kwdef struct ViscousSponge{FT} <: AbstractSponge
@@ -560,6 +562,7 @@ Base.@kwdef struct AtmosModel{
     rayleigh_sponge::RS = nothing
     sfc_temperature::ST = nothing
     insolation::IN = nothing
+    disable_surface_flux_tendency::Bool = false
     surface_model::SM = nothing
     surface_albedo::SA = nothing
     numerics::NUM = nothing
diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl
index cf92e39292..b6687704f8 100644
--- a/src/surface_conditions/surface_setups.jl
+++ b/src/surface_conditions/surface_setups.jl
@@ -1,5 +1,13 @@
 abstract type SurfaceSetup end
 
+"""
+    NoSurface()
+
+Used to indicate that there is no surface flux.
+"""
+struct NoSurface <: SurfaceSetup end
+(surface_setup::NoSurface)(params) = nothing
+
 """
     PrescribedSurface()