diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 98bca456ec7..50a5a6f62e7 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -723,6 +723,58 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml
         artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist_edmfx/*"
+        soft_fail: true
+        agents:
+          slurm_mem: 20GB
+      
+      - label: ":genie: Advective EDMFX advection test in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_adv_test_box.yml
+        artifact_paths: "advective_edmfx_adv_test_box/*"
+        agents:
+          slurm_mem: 20GB
+
+      - label: ":genie: Advective EDMFX GABLS in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_gabls_box.yml
+        artifact_paths: "advective_edmfx_gabls_box/*"
+        agents:
+          slurm_mem: 20GB
+
+      - label: ":genie: Advective EDMFX Bomex in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_bomex_box.yml
+        artifact_paths: "advective_edmfx_bomex_box/*"
+        soft_fail: true
+        agents:
+          slurm_mem: 20GB
+
+      - label: ":genie: Advective EDMFX Dycoms RF01 in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_dycoms_rf01_box.yml
+        artifact_paths: "advective_edmfx_dycoms_rf01_box/*"
+        agents:
+          slurm_mem: 20GB
+      
+      - label: ":genie: Advective EDMFX Rico in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_rico_box.yml
+        artifact_paths: "advective_edmfx_rico_box/*"
+        soft_fail: true
+        agents:
+          slurm_mem: 20GB
+
+      - label: ":genie: Advective EDMFX TRMM in a box"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/advective_edmfx_trmm_box.yml
+        artifact_paths: "advective_edmfx_trmm_box/*"
+        soft_fail: true
         agents:
           slurm_mem: 20GB
 
diff --git a/config/longrun_configs/longrun_aquaplanet_amip.yml b/config/longrun_configs/longrun_aquaplanet_amip.yml
index 4251539b34a..b7c7d20b9a2 100644
--- a/config/longrun_configs/longrun_aquaplanet_amip.yml
+++ b/config/longrun_configs/longrun_aquaplanet_amip.yml
@@ -14,7 +14,7 @@ turbconv: "diagnostic_edmfx"
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/advective_edmfx_adv_test_box.yml b/config/model_configs/advective_edmfx_adv_test_box.yml
new file mode 100644
index 00000000000..637640b61c5
--- /dev/null
+++ b/config/model_configs/advective_edmfx_adv_test_box.yml
@@ -0,0 +1,23 @@
+job_id: "advective_edmfx_adv_test_box"
+initial_condition: "MoistAdiabaticProfileEDMFX"
+advection_test: true
+turbconv: "advective_edmfx"
+ode_algo: "SSP33ShuOsher"
+edmfx_upwinding: "first_order"
+config: "box"
+moist: "equil"
+hyperdiff: "true"
+kappa_4: 1e8
+x_max: 1e4
+y_max: 1e4
+z_max: 5.5e4
+x_elem: 2
+y_elem: 2
+z_elem: 63
+dz_bottom: 30.0
+dz_top: 3000.0
+dt: "10secs"
+t_end: "3600secs"
+dt_save_to_disk: "100secs"
+FLOAT_TYPE: "Float64"
+toml: [toml/edmfx_box_advection.toml]
diff --git a/config/model_configs/advective_edmfx_bomex_box.yml b/config/model_configs/advective_edmfx_bomex_box.yml
new file mode 100644
index 00000000000..9c611f900e3
--- /dev/null
+++ b/config/model_configs/advective_edmfx_bomex_box.yml
@@ -0,0 +1,30 @@
+job_id: "advective_edmfx_bomex_box"
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
+surface_setup: "Bomex"
+turbconv: "advective_edmfx"
+edmfx_upwinding: first_order
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+prognostic_tke: false
+moist: "equil"
+config: "box"
+hyperdiff: "true"
+kappa_4: 1.0e12
+x_max: 1e5
+y_max: 1e5
+z_max: 3e3
+x_elem: 2
+y_elem: 2
+z_elem: 60
+z_stretch: false
+perturb_initstate: false
+dt: "5secs"
+t_end: "6hours"
+dt_save_to_disk: "5mins"
+toml: [toml/advective_edmfx_box.toml]
diff --git a/config/model_configs/advective_edmfx_dycoms_rf01_box.yml b/config/model_configs/advective_edmfx_dycoms_rf01_box.yml
new file mode 100644
index 00000000000..3c59eab1b73
--- /dev/null
+++ b/config/model_configs/advective_edmfx_dycoms_rf01_box.yml
@@ -0,0 +1,29 @@
+job_id: advective_edmfx_dycoms_rf01_box
+initial_condition: DYCOMS_RF01
+subsidence: DYCOMS
+edmf_coriolis: DYCOMS_RF01
+rad: DYCOMS_RF01
+surface_setup: DYCOMS_RF01
+turbconv: advective_edmfx
+edmfx_upwinding: first_order
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+prognostic_tke: false
+moist: equil
+config: box
+hyperdiff: "true"
+kappa_4: 1e12
+x_max: 1e5
+y_max: 1e5
+z_max: 1500
+x_elem: 2
+y_elem: 2
+z_elem: 30
+z_stretch: false
+dt: 5secs
+t_end: 4hours
+dt_save_to_disk: 10mins
+toml: [toml/advective_edmfx_box.toml]
diff --git a/config/model_configs/advective_edmfx_gabls_box.yml b/config/model_configs/advective_edmfx_gabls_box.yml
new file mode 100644
index 00000000000..50ba179be36
--- /dev/null
+++ b/config/model_configs/advective_edmfx_gabls_box.yml
@@ -0,0 +1,29 @@
+job_id: "advective_edmfx_gabls_box"
+initial_condition: GABLS
+edmf_coriolis: GABLS
+surface_setup: GABLS
+turbconv: "advective_edmfx"
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+prognostic_tke: false
+moist: "equil"
+config: "box"
+hyperdiff: "true"
+kappa_4: 1e12
+x_max: 1e5
+y_max: 1e5
+z_max: 400
+x_elem: 2
+y_elem: 2
+z_elem: 8
+z_stretch: false
+dt: "5secs"
+t_end: "9hours"
+dt_save_to_disk: "10mins"
+perturb_initstate: false
+FLOAT_TYPE: "Float64"
+toml: [toml/advective_edmfx_box.toml]
diff --git a/config/model_configs/advective_edmfx_rico_box.yml b/config/model_configs/advective_edmfx_rico_box.yml
new file mode 100644
index 00000000000..0320db52d38
--- /dev/null
+++ b/config/model_configs/advective_edmfx_rico_box.yml
@@ -0,0 +1,30 @@
+job_id: "advective_edmfx_rico_box"
+initial_condition: "Rico"
+subsidence: "Rico"
+ls_adv: "Rico"
+surface_setup: "Rico"
+turbconv: "advective_edmfx"
+edmfx_upwinding: first_order
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+prognostic_tke: false
+moist: "equil"
+precip_model: "0M"
+config: "box"
+hyperdiff: "true"
+kappa_4: 1e12
+x_max: 1e5
+y_max: 1e5
+z_max: 4e3
+x_elem: 2
+y_elem: 2
+z_elem: 100
+z_stretch: false
+perturb_initstate: false
+dt: "5secs"
+t_end: "6hours"
+dt_save_to_disk: "10mins"
+toml: [toml/advective_edmfx_box.toml]
diff --git a/config/model_configs/advective_edmfx_trmm_box.yml b/config/model_configs/advective_edmfx_trmm_box.yml
new file mode 100644
index 00000000000..fc062af52d9
--- /dev/null
+++ b/config/model_configs/advective_edmfx_trmm_box.yml
@@ -0,0 +1,30 @@
+job_id: advective_edmfx_trmm_box
+initial_condition: TRMM_LBA
+rad: TRMM_LBA
+surface_setup: TRMM_LBA
+turbconv: advective_edmfx
+edmfx_upwinding: first_order
+edmfx_entr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+prognostic_tke: false
+moist: equil
+apply_limiter: false
+precip_model: "0M"
+config: box
+hyperdiff: "true"
+kappa_4: 1e12
+x_max: 1e5
+y_max: 1e5
+z_max: 16400
+x_elem: 2
+y_elem: 2
+z_elem: 82
+z_stretch: false
+dt: 5secs
+t_end: 6hours
+dt_save_to_disk: 10mins
+FLOAT_TYPE: "Float64"
+toml: [toml/advective_edmfx_box.toml]
diff --git a/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml b/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml
index c1c2a37bf73..a1fdf42e5d3 100644
--- a/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml
+++ b/config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml
@@ -6,7 +6,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order 
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_bomex_box.yml b/config/model_configs/diagnostic_edmfx_bomex_box.yml
index 6f172eb58e4..0aed107f8f6 100644
--- a/config/model_configs/diagnostic_edmfx_bomex_box.yml
+++ b/config/model_configs/diagnostic_edmfx_bomex_box.yml
@@ -9,7 +9,7 @@ turbconv: "diagnostic_edmfx"
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
index 23ebd68de84..9a50a12adef 100644
--- a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
+++ b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
@@ -9,7 +9,7 @@ turbconv: "diagnostic_edmfx"
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
index de7af5b0f81..a7a26db7c22 100644
--- a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
+++ b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
@@ -8,7 +8,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order 
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_gabls_box.yml b/config/model_configs/diagnostic_edmfx_gabls_box.yml
index 74c38a0b3af..4feb77db8ce 100644
--- a/config/model_configs/diagnostic_edmfx_gabls_box.yml
+++ b/config/model_configs/diagnostic_edmfx_gabls_box.yml
@@ -6,7 +6,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true 
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_rico_box.yml b/config/model_configs/diagnostic_edmfx_rico_box.yml
index ef7d30f2a8c..77484647164 100644
--- a/config/model_configs/diagnostic_edmfx_rico_box.yml
+++ b/config/model_configs/diagnostic_edmfx_rico_box.yml
@@ -8,7 +8,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml
index 2ccbc4da6c9..d35e23799d7 100644
--- a/config/model_configs/diagnostic_edmfx_trmm_box.yml
+++ b/config/model_configs/diagnostic_edmfx_trmm_box.yml
@@ -7,7 +7,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
index 32f585b6ccb..5bbd68dc934 100644
--- a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
+++ b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
@@ -7,7 +7,7 @@ turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantTimescale" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
diff --git a/config/model_configs/edmfx_bomex_box.yml b/config/model_configs/edmfx_bomex_box.yml
index 8aed4c6a7c3..507a8528c39 100644
--- a/config/model_configs/edmfx_bomex_box.yml
+++ b/config/model_configs/edmfx_bomex_box.yml
@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
diff --git a/config/model_configs/edmfx_dycoms_rf01_box.yml b/config/model_configs/edmfx_dycoms_rf01_box.yml
index c59a9f7158b..1da607d5ac0 100644
--- a/config/model_configs/edmfx_dycoms_rf01_box.yml
+++ b/config/model_configs/edmfx_dycoms_rf01_box.yml
@@ -7,7 +7,7 @@ surface_setup: DYCOMS_RF01
 turbconv: edmfx
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -23,7 +23,7 @@ x_elem: 2
 y_elem: 2
 z_elem: 30
 z_stretch: false
-dt: 2secs
+dt: 5secs
 t_end: 4hours
-dt_save_to_disk: 10secs
+dt_save_to_disk: 10mins
 toml: [toml/edmfx_box.toml]
diff --git a/config/model_configs/edmfx_gabls_box.yml b/config/model_configs/edmfx_gabls_box.yml
index c61d11f7b14..b01eec483b2 100644
--- a/config/model_configs/edmfx_gabls_box.yml
+++ b/config/model_configs/edmfx_gabls_box.yml
@@ -5,7 +5,7 @@ surface_setup: GABLS
 turbconv: "edmfx"
 edmfx_upwinding: "first_order"
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
diff --git a/config/model_configs/edmfx_rico_box.yml b/config/model_configs/edmfx_rico_box.yml
index af02f49bb3b..634a1a91d3b 100644
--- a/config/model_configs/edmfx_rico_box.yml
+++ b/config/model_configs/edmfx_rico_box.yml
@@ -6,7 +6,7 @@ surface_setup: "Rico"
 turbconv: "edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
diff --git a/config/model_configs/edmfx_trmm_box.yml b/config/model_configs/edmfx_trmm_box.yml
index 3720fec2f38..2b790f5f316 100644
--- a/config/model_configs/edmfx_trmm_box.yml
+++ b/config/model_configs/edmfx_trmm_box.yml
@@ -5,7 +5,7 @@ surface_setup: TRMM_LBA
 turbconv: edmfx
 edmfx_upwinding: first_order
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
diff --git a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml b/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml
index 656cb478f72..0d2d839d735 100644
--- a/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml
+++ b/config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml
@@ -2,7 +2,7 @@ dt_save_to_disk: "30secs"
 initial_condition: "MoistBaroclinicWaveWithEDMF"
 dt: "1secs"
 edmfx_entr_model: "ConstantCoefficient" 
-edmfx_detr_model: "ConstantCoefficient" 
+edmfx_detr_model: "BOverW" 
 t_end: "6mins"
 turbconv: "edmfx"
 ode_algo: "SSP33ShuOsher"
diff --git a/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml b/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml
index b5f20d92f3b..e67d1b10b4a 100644
--- a/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml
+++ b/config/perf_configs/flame_perf_target_diagnostic_edmfx.yml
@@ -6,4 +6,4 @@ prognostic_tke: true
 edmfx_sgs_flux: true
 turbconv: "diagnostic_edmfx"
 job_id: "flame_perf_target_diagnostic_edmfx"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
diff --git a/config/perf_configs/flame_perf_target_edmfx.yml b/config/perf_configs/flame_perf_target_edmfx.yml
index f892c20bd42..0d5a514bc44 100644
--- a/config/perf_configs/flame_perf_target_edmfx.yml
+++ b/config/perf_configs/flame_perf_target_edmfx.yml
@@ -3,7 +3,7 @@ reference_job_id: "sphere_baroclinic_wave_rhoe_equilmoist"
 ode_algo: "SSP33ShuOsher"
 initial_condition: "MoistBaroclinicWaveWithEDMF"
 edmfx_entr_model: "ConstantCoefficient"
-edmfx_detr_model: "ConstantCoefficient"
+edmfx_detr_model: "BOverW"
 edmfx_sgs_flux: true
 turbconv: "edmfx"
 job_id: "flame_perf_target_edmfx"
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 6af89053b03..552c94041ad 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -40,6 +40,7 @@ reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id
 
 is_edmfx =
     atmos.turbconv_model isa CA.EDMFX ||
+    atmos.turbconv_model isa CA.AdvectiveEDMFX ||
     atmos.turbconv_model isa CA.DiagnosticEDMFX
 if is_edmfx && config.parsed_args["post_process"]
     contours_and_profiles(simulation.output_dir, reference_job_id)
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
index 1f0c7fcac38..3daa52df33c 100644
--- a/src/ClimaAtmos.jl
+++ b/src/ClimaAtmos.jl
@@ -26,6 +26,7 @@ include(joinpath("TurbulenceConvection_deprecated", "TurbulenceConvection.jl"))
 import .TurbulenceConvection as TC
 
 include(joinpath("cache", "edmf_precomputed_quantities.jl"))
+include(joinpath("cache", "advective_edmf_precomputed_quantities.jl"))
 include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl"))
 include(joinpath("cache", "precomputed_quantities.jl"))
 
diff --git a/src/cache/advective_edmf_precomputed_quantities.jl b/src/cache/advective_edmf_precomputed_quantities.jl
new file mode 100644
index 00000000000..b47f81e0359
--- /dev/null
+++ b/src/cache/advective_edmf_precomputed_quantities.jl
@@ -0,0 +1,268 @@
+#####
+##### Precomputed quantities
+#####
+import Thermodynamics as TD
+import ClimaCore: Spaces, Fields
+
+"""
+    set_advective_edmf_precomputed_quantities!(Y, p, t)
+
+Updates the precomputed quantities stored in `p` for edmfx.
+"""
+function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
+    (; energy_form, moisture_model, turbconv_model) = p.atmos
+    #EDMFX BCs only support total energy as state variable
+    @assert energy_form isa TotalEnergy
+    @assert !(moisture_model isa DryModel)
+
+    FT = Spaces.undertype(axes(Y.c))
+    (; params) = p
+    (; dt) = p.simulation
+    thermo_params = CAP.thermodynamics_params(params)
+    n = n_mass_flux_subdomains(turbconv_model)
+    thermo_args = (thermo_params, energy_form, moisture_model)
+
+    (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p
+    (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
+    (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p
+    (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
+
+    @. ᶜρa⁰ = ρa⁰(Y.c)
+    @. ᶜtke⁰ = divide_by_ρa(Y.c.sgs⁰.ρatke, ᶜρa⁰, 0, Y.c.ρ, turbconv_model)
+    # TODO: use all updrafts
+    @. ᶜh_tot⁰ = divide_by_ρa(
+        Y.c.ρ * ᶜh_tot - (Y.c.sgsʲs.:1).ρa * (Y.c.sgsʲs.:1).h_tot,
+        ᶜρa⁰,
+        Y.c.ρ * ᶜh_tot,
+        Y.c.ρ,
+        turbconv_model,
+    )
+    @. ᶜq_tot⁰ = divide_by_ρa(
+        Y.c.ρq_tot - (Y.c.sgsʲs.:1).ρa * (Y.c.sgsʲs.:1).q_tot,
+        ᶜρa⁰,
+        Y.c.ρq_tot,
+        Y.c.ρ,
+        turbconv_model,
+    )
+    set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
+    set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
+    @. ᶜK⁰ += ᶜtke⁰
+    @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ, ᶜq_tot⁰)
+    @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
+
+    for j in 1:n
+        ᶜuʲ = ᶜuʲs.:($j)
+        ᶠu³ʲ = ᶠu³ʲs.:($j)
+        ᶜKʲ = ᶜKʲs.:($j)
+        ᶠu₃ʲ = Y.f.sgsʲs.:($j).u₃
+        ᶜtsʲ = ᶜtsʲs.:($j)
+        ᶜρʲ = ᶜρʲs.:($j)
+        ᶜh_totʲ = Y.c.sgsʲs.:($j).h_tot
+        ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
+
+        set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
+        @. ᶜtsʲ =
+            TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_totʲ - ᶜKʲ - ᶜΦ, ᶜq_totʲ)
+        @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
+
+        # When ᶜe_intʲ = ᶜe_int and ᶜq_totʲ = ᶜq_tot, we still observe that
+        # ᶜρʲ != ᶜρ. This is because the conversion from ᶜρ to ᶜp to ᶜρʲ
+        # introduces a tiny round-off error of order epsilon to ᶜρʲ. If left
+        # unchecked, this round-off error then changes the tendency of ᶠu₃ʲ,
+        # which in turn introduces an error to ᶠu³ʲ, which then increases
+        # the error in ᶜρʲ. For now, we will filter ᶜρʲ to fix this. Note
+        # that this will no longer be necessary after we add diffusion.
+        # @. ᶜρʲ = ifelse(abs(ᶜρʲ - Y.c.ρ) <= 2 * eps(Y.c.ρ), Y.c.ρ, ᶜρʲ)
+
+        # EDMFX boundary condition:
+
+        # We need field_values everywhere because we are mixing
+        # information from surface and first interior inside the
+        # sgs_h/q_tot_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.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 ᶜh_totʲ...
+        ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
+        ᶜh_totʲ_int_val = Fields.field_values(Fields.level(ᶜh_totʲ, 1))
+        @. ᶜh_totʲ_int_val = sgs_scalar_first_interior_bc(
+            ᶜz_int_val - z_sfc_val,
+            ᶜρ_int_val,
+            ᶜh_tot_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_int_val = Fields.field_values(Fields.level(ᶜspecific.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,
+            ᶜ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.
+        ᶜKʲ_int_val = Fields.field_values(Fields.level(ᶜKʲ, 1))
+        ᶜΦ_int_val = Fields.field_values(Fields.level(ᶜΦ, 1))
+        ᶜtsʲ_int_val = Fields.field_values(Fields.level(ᶜtsʲ, 1))
+        @. ᶜtsʲ_int_val = TD.PhaseEquil_phq(
+            thermo_params,
+            ᶜp_int_val,
+            ᶜh_totʲ_int_val - ᶜKʲ_int_val - ᶜΦ_int_val,
+            ᶜq_totʲ_int_val,
+        )
+        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))
+
+        turbconv_params = CAP.turbconv_params(params)
+        @. sgsʲs_ρa_int_val =
+            $(FT(turbconv_params.surface_area)) *
+            TD.air_density(thermo_params, ᶜtsʲ_int_val)
+    end
+
+    ᶜz = Fields.coordinate_field(Y.c).z
+    z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+    ᶜdz = Fields.Δz_field(axes(Y.c))
+    ᶜlg = Fields.local_geometry_field(Y.c)
+
+    for j in 1:n
+        @. ᶜentrʲs.:($$j) = entrainment(
+            params,
+            ᶜz,
+            z_sfc,
+            ᶜp,
+            Y.c.ρ,
+            buoyancy_flux,
+            draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
+            get_physical_w(ᶜuʲs.:($$j), ᶜlg),
+            TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
+            ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρʲs.:($$j)),
+            get_physical_w(ᶜu⁰, ᶜlg),
+            TD.relative_humidity(thermo_params, ᶜts⁰),
+            ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰),
+            dt,
+            p.atmos.edmfx_entr_model,
+        )
+        @. ᶜdetrʲs.:($$j) = detrainment(
+            params,
+            ᶜz,
+            z_sfc,
+            ᶜp,
+            Y.c.ρ,
+            buoyancy_flux,
+            draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
+            get_physical_w(ᶜuʲs.:($$j), ᶜlg),
+            TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
+            ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρʲs.:($$j)),
+            get_physical_w(ᶜu⁰, ᶜlg),
+            TD.relative_humidity(thermo_params, ᶜts⁰),
+            ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰),
+            dt,
+            p.atmos.edmfx_detr_model,
+        )
+    end
+
+    # First order approximation: Use environmental mean fields.
+    @. ᶜlinear_buoygrad = buoyancy_gradients(
+        params,
+        moisture_model,
+        EnvBuoyGrad(
+            BuoyGradMean(),
+            TD.air_temperature(thermo_params, ᶜts⁰),                           # t_sat
+            TD.vapor_specific_humidity(thermo_params, ᶜts⁰),                   # qv_sat
+            ᶜq_tot⁰,                                                           # qt_sat
+            TD.dry_pottemp(thermo_params, ᶜts⁰),                               # θ_sat
+            TD.liquid_ice_pottemp(thermo_params, ᶜts⁰),                        # θ_liq_ice_sat
+            projected_vector_data(
+                C3,
+                ᶜgradᵥ(ᶠinterp(TD.virtual_pottemp(thermo_params, ᶜts⁰))),
+                ᶜlg,
+            ),                                                                 # ∂θv∂z_unsat
+            projected_vector_data(C3, ᶜgradᵥ(ᶠinterp(ᶜq_tot⁰)), ᶜlg),          # ∂qt∂z_sat
+            projected_vector_data(
+                C3,
+                ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts⁰))),
+                ᶜlg,
+            ),                                                                 # ∂θl∂z_sat
+            ᶜp,                                                                # p
+            ifelse(TD.has_condensate(thermo_params, ᶜts⁰), 1, 0),              # en_cld_frac
+            ᶜρ⁰,                                                               # ρ
+        ),
+    )
+
+    # TODO: Currently the shear production only includes vertical gradients
+    ᶠu⁰ = p.ᶠtemp_C123
+    @. ᶠu⁰ = C123(ᶠinterp(Y.c.uₕ)) + C123(ᶠu³⁰)
+    ᶜstrain_rate = p.ᶜtemp_UVWxUVW
+    compute_strain_rate_center!(ᶜstrain_rate, ᶠu⁰)
+    @. ᶜstrain_rate_norm = norm_sqr(ᶜstrain_rate)
+
+    ᶜprandtl_nvec = p.ᶜtemp_scalar
+    @. ᶜprandtl_nvec = turbulent_prandtl_number(
+        params,
+        obukhov_length,
+        ᶜlinear_buoygrad,
+        ᶜstrain_rate_norm,
+    )
+    ᶜtke_exch = p.ᶜtemp_scalar_2
+    @. ᶜtke_exch = 0
+    for j in 1:n
+        @. ᶜtke_exch +=
+            Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * (
+                1 / 2 *
+                (
+                    get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg)
+                )^2 - ᶜtke⁰
+            )
+    end
+
+    sfc_tke = Fields.level(ᶜtke⁰, 1)
+    @. ᶜmixing_length = mixing_length(
+        p.params,
+        ustar,
+        ᶜz,
+        z_sfc,
+        ᶜdz,
+        sfc_tke,
+        ᶜlinear_buoygrad,
+        ᶜtke⁰,
+        obukhov_length,
+        ᶜstrain_rate_norm,
+        ᶜprandtl_nvec,
+        ᶜtke_exch,
+    )
+
+    turbconv_params = CAP.turbconv_params(params)
+    c_m = TCP.tke_ed_coeff(turbconv_params)
+    @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0))
+    # TODO: add Prantdl number
+    @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec
+
+    return nothing
+end
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index b3be3719ead..0349a4b0896 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -40,6 +40,10 @@ function precomputed_quantities(Y, atmos)
         !(atmos.moisture_model isa DryModel) &&
         atmos.energy_form isa TotalEnergy
     ) || !(atmos.turbconv_model isa DiagnosticEDMFX)
+    @assert (
+        !(atmos.moisture_model isa DryModel) &&
+        atmos.energy_form isa TotalEnergy
+    ) || !(atmos.turbconv_model isa AdvectiveEDMFX)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     n = n_mass_flux_subdomains(atmos.turbconv_model)
@@ -88,6 +92,32 @@ function precomputed_quantities(Y, atmos)
                 ) : (;)
             )...,
         ) : (;)
+    advective_sgs_quantities =
+        atmos.turbconv_model isa AdvectiveEDMFX ?
+        (;
+            ᶜtke⁰ = similar(Y.c, FT),
+            ᶜρa⁰ = similar(Y.c, FT),
+            ᶠu₃⁰ = similar(Y.f, C3{FT}),
+            ᶜu⁰ = similar(Y.c, C123{FT}),
+            ᶠu³⁰ = similar(Y.f, CT3{FT}),
+            ᶜK⁰ = similar(Y.c, FT),
+            ᶜh_tot⁰ = similar(Y.c, FT),
+            ᶜq_tot⁰ = similar(Y.c, FT),
+            ᶜts⁰ = similar(Y.c, TST),
+            ᶜρ⁰ = similar(Y.c, FT),
+            ᶜlinear_buoygrad = similar(Y.c, FT),
+            ᶜstrain_rate_norm = similar(Y.c, FT),
+            ᶜmixing_length = similar(Y.c, FT),
+            ᶜK_u = similar(Y.c, FT),
+            ᶜK_h = similar(Y.c, FT),
+            ᶜuʲs = similar(Y.c, NTuple{n, C123{FT}}),
+            ᶠu³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
+            ᶜKʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
+            ᶜρʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
+        ) : (;)
     diagnostic_sgs_quantities =
         atmos.turbconv_model isa DiagnosticEDMFX ?
         (;
@@ -117,7 +147,12 @@ function precomputed_quantities(Y, atmos)
             ᶜK_h = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
         ) : (;)
-    return (; gs_quantities..., sgs_quantities..., diagnostic_sgs_quantities...)
+    return (;
+        gs_quantities...,
+        sgs_quantities...,
+        advective_sgs_quantities...,
+        diagnostic_sgs_quantities...,
+    )
 end
 
 # Interpolates the third contravariant component of Y.c.uₕ to cell faces.
@@ -140,7 +175,7 @@ function set_velocity_at_surface!(Y, ᶠuₕ³, turbconv_model)
     sfc_uₕ³ = Fields.level(ᶠuₕ³.components.data.:1, half)
     sfc_g³³ = g³³_field(sfc_u₃)
     @. sfc_u₃ = -sfc_uₕ³ / sfc_g³³ # u³ = uₕ³ + w³ = uₕ³ + w₃ * g³³
-    if turbconv_model isa EDMFX
+    if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX
         for j in 1:n_mass_flux_subdomains(turbconv_model)
             sfc_u₃ʲ = Fields.level(Y.f.sgsʲs.:($j).u₃.components.data.:1, half)
             @. sfc_u₃ʲ = sfc_u₃
@@ -160,7 +195,7 @@ function set_velocity_at_top!(Y, turbconv_model)
         Spaces.nlevels(axes(Y.c)) + half,
     )
     @. top_u₃ = 0
-    if turbconv_model isa EDMFX
+    if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX
         for j in 1:n_mass_flux_subdomains(turbconv_model)
             top_u₃ʲ = Fields.level(
                 Y.f.sgsʲs.:($j).u₃.components.data.:1,
@@ -341,6 +376,10 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
         set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     end
 
+    if turbconv_model isa AdvectiveEDMFX
+        set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
+    end
+
     if turbconv_model isa DiagnosticEDMFX
         set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
     end
@@ -373,6 +412,27 @@ function output_sgs_quantities(Y, p, t)
     return (; ᶜspecific⁺, ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰)
 end
 
+"""
+    output_advective_sgs_quantities(Y, p, t)
+
+Sets `ᶜu⁺`, `ᶠu³⁺`, `ᶜts⁺` and `ᶜa⁺` to be the same as the
+values of the first updraft.
+"""
+function output_advective_sgs_quantities(Y, p, t)
+    (; turbconv_model) = p.atmos
+    thermo_params = CAP.thermodynamics_params(p.params)
+    (; ᶜp, ᶜρa⁰, ᶜρ⁰, ᶜΦ, ᶜtsʲs) = p
+    ᶠuₕ³ = p.ᶠtemp_CT3
+    set_ᶠuₕ³!(ᶠuₕ³, Y)
+    (ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺) = similar.((p.ᶠu₃⁰, p.ᶜu⁰, p.ᶠu³⁰, p.ᶜK⁰))
+    set_sgs_ᶠu₃!(u₃⁺, ᶠu₃⁺, Y, turbconv_model)
+    set_velocity_quantities!(ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶠu₃⁺, Y.c.uₕ, ᶠuₕ³)
+    ᶜts⁺ = ᶜtsʲs.:1
+    ᶜa⁺ = @. draft_area(ρa⁺(Y.c), TD.air_density(thermo_params, ᶜts⁺))
+    ᶜa⁰ = @. draft_area(ᶜρa⁰, ᶜρ⁰)
+    return (; ᶠu₃⁺, ᶜu⁺, ᶠu³⁺, ᶜK⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰)
+end
+
 """
     output_diagnostic_sgs_quantities(Y, p, t)
 
diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl
index 4625e50872c..4112e8d7e76 100644
--- a/src/callbacks/callbacks.jl
+++ b/src/callbacks/callbacks.jl
@@ -304,6 +304,37 @@ NVTX.@annotate function compute_diagnostics(integrator)
                 ᶜts⁰,
             ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
         )
+    elseif p.atmos.turbconv_model isa AdvectiveEDMFX
+        (; ᶜtke⁰, ᶜu⁰, ᶜts⁰, ᶜmixing_length) = p
+        (; ᶜu⁺, ᶜts⁺, ᶜa⁺, ᶜa⁰) = output_advective_sgs_quantities(Y, p, t)
+        env_diagnostics = (;
+            common_diagnostics(p, ᶜu⁰, ᶜts⁰)...,
+            area = ᶜa⁰,
+            cloud_fraction = get_cloud_fraction.(
+                thermo_params,
+                env_thermo_quad,
+                ᶜp,
+                ᶜts⁰,
+            ),
+            tke = ᶜtke⁰,
+            mixing_length = ᶜmixing_length,
+        )
+        draft_diagnostics = (;
+            common_diagnostics(p, ᶜu⁺, ᶜts⁺)...,
+            area = ᶜa⁺,
+            cloud_fraction = cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
+        )
+        turbulence_convection_diagnostic = (;
+            add_prefix(env_diagnostics, :env_)...,
+            add_prefix(draft_diagnostics, :draft_)...,
+            cloud_fraction = ᶜa⁰ .*
+                             get_cloud_fraction.(
+                thermo_params,
+                env_thermo_quad,
+                ᶜp,
+                ᶜts⁰,
+            ) .+ ᶜa⁺ .* cloud_fraction.(thermo_params, ᶜts⁺, ᶜa⁺),
+        )
     elseif p.atmos.turbconv_model isa DiagnosticEDMFX
         (; ᶜtke⁰, ᶜmixing_length) = p
         (; ᶜu⁺, ᶜts⁺, ᶜa⁺) = output_diagnostic_sgs_quantities(Y, p, t)
diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl
index 29cef59fd7c..081bc892a18 100644
--- a/src/initial_conditions/InitialConditions.jl
+++ b/src/initial_conditions/InitialConditions.jl
@@ -16,6 +16,7 @@ import ..PrognosticSurfaceTemperature
 import ..C3
 import ..C12
 import ..EDMFX
+import ..AdvectiveEDMFX
 import ..DiagnosticEDMFX
 import ..n_mass_flux_subdomains
 import ..gs_to_sgs
diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl
index 84865dd0d24..8caf9e21705 100644
--- a/src/initial_conditions/atmos_state.jl
+++ b/src/initial_conditions/atmos_state.jl
@@ -129,6 +129,7 @@ function turbconv_center_variables(ls, turbconv_model::TC.EDMFModel, gs_vars)
     up = ntuple(_ -> (; ρarea = ρa, ρae_tot, ρaq_tot), Val(n))
     return (; turbconv = (; en, up))
 end
+
 function turbconv_center_variables(ls, turbconv_model::EDMFX, gs_vars)
     n = n_mass_flux_subdomains(turbconv_model)
     a_draft = ls.turbconv_state.draft_area
@@ -137,6 +138,20 @@ function turbconv_center_variables(ls, turbconv_model::EDMFX, gs_vars)
     return (; sgs⁰, sgsʲs)
 end
 
+function turbconv_center_variables(ls, turbconv_model::AdvectiveEDMFX, gs_vars)
+    n = n_mass_flux_subdomains(turbconv_model)
+    a_draft = ls.turbconv_state.draft_area
+    sgs⁰ = (; ρatke = ls.ρ * (1 - a_draft) * ls.turbconv_state.tke)
+    ρa = ls.ρ * a_draft / n
+    h_tot =
+        TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) +
+        norm_sqr(ls.velocity) / 2 +
+        CAP.grav(ls.params) * ls.geometry.coordinates.z
+    q_tot = TD.total_specific_humidity(ls.thermo_params, ls.thermo_state)
+    sgsʲs = ntuple(_ -> (; ρa = ρa, h_tot = h_tot, q_tot = q_tot), Val(n))
+    return (; sgs⁰, sgsʲs)
+end
+
 function turbconv_center_variables(ls, turbconv_model::DiagnosticEDMFX, gs_vars)
     sgs⁰ = (; ρatke = ls.ρ * ls.turbconv_state.tke)
     return (; sgs⁰)
@@ -157,4 +172,10 @@ turbconv_face_variables(ls, turbconv_model::EDMFX) = (;
         Val(n_mass_flux_subdomains(turbconv_model)),
     )
 )
+turbconv_face_variables(ls, turbconv_model::AdvectiveEDMFX) = (;
+    sgsʲs = ntuple(
+        _ -> (; u₃ = C3(ls.turbconv_state.velocity, ls.geometry)),
+        Val(n_mass_flux_subdomains(turbconv_model)),
+    )
+)
 turbconv_face_variables(ls, turbconv_model::DiagnosticEDMFX) = (;)
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index 94d4bbfb496..84e336ed010 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -12,12 +12,14 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
     if p.atmos.turbconv_model isa AbstractEDMF
         (; ᶜu⁰) = p
     end
-    if p.atmos.turbconv_model isa EDMFX
+    if p.atmos.turbconv_model isa EDMFX ||
+       p.atmos.turbconv_model isa AdvectiveEDMFX
         (; ᶜuʲs) = p
     end
 
     @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu)
-    if p.atmos.turbconv_model isa EDMFX
+    if p.atmos.turbconv_model isa EDMFX ||
+       p.atmos.turbconv_model isa AdvectiveEDMFX
         for j in 1:n
             @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j))
         end
@@ -42,6 +44,15 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
         end
     end
 
+    if p.atmos.turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            @. Yₜ.c.sgsʲs.:($$j).h_tot -=
+                adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).h_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).q_tot)
+        end
+    end
+
     if use_prognostic_tke(p.atmos.turbconv_model)
         @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰)
     end
@@ -71,6 +82,14 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
             end
         end
     end
+
+    # FIXME: add horizontal advection
+    if p.atmos.turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -= 0
+        end
+    end
+
     return nothing
 end
 
@@ -84,13 +103,16 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     ᶜJ = Fields.local_geometry_field(Y.c).J
     (; ᶜu, ᶠu³, ᶜK, ᶜf) = p
     (; edmfx_upwinding) = n > 0 || advect_tke ? p : all_nothing
-    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs, ᶜspecificʲs) = n > 0 ? p : all_nothing
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = n > 0 ? p : all_nothing
+    (; ᶜspecificʲs) = turbconv_model isa EDMFX ? p : all_nothing
     (; ᶜp, ᶜp_ref, ᶜρ_ref, ᶠgradᵥ_ᶜΦ) = n > 0 ? p : all_nothing
-    (; ᶜh_totʲs) = n > 0 && is_total_energy ? p : all_nothing
+    (; ᶜh_totʲs) = turbconv_model isa EDMFX && is_total_energy ? p : all_nothing
     (; ᶠu³⁰) = advect_tke ? p : all_nothing
     ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing
     ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing
-    ᶜtke⁰ = advect_tke ? (n > 0 ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing
+    ᶜtke⁰ =
+        advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) :
+        nothing
     ᶜa_scalar = p.ᶜtemp_scalar
     ᶜω³ = p.ᶜtemp_CT3
     ᶠω¹² = p.ᶠtemp_CT12
@@ -139,25 +161,12 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         end
 
         # TODO: Move this to implicit_vertical_advection_tendency!.
-        for j in 1:n
-            @. ᶜa_scalar[colidx] =
-                draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx])
-            vertical_transport!(
-                Yₜ.c.sgsʲs.:($j).ρa[colidx],
-                ᶜJ[colidx],
-                ᶜρʲs.:($j)[colidx],
-                ᶠu³ʲs.:($j)[colidx],
-                ᶜa_scalar[colidx],
-                dt,
-                edmfx_upwinding,
-            )
-
-            if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j))
+        if p.atmos.turbconv_model isa EDMFX
+            for j in 1:n
                 @. ᶜa_scalar[colidx] =
-                    ᶜh_totʲs.:($$j)[colidx] *
                     draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx])
                 vertical_transport!(
-                    Yₜ.c.sgsʲs.:($j).ρae_tot[colidx],
+                    Yₜ.c.sgsʲs.:($j).ρa[colidx],
                     ᶜJ[colidx],
                     ᶜρʲs.:($j)[colidx],
                     ᶠu³ʲs.:($j)[colidx],
@@ -165,16 +174,51 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     dt,
                     edmfx_upwinding,
                 )
+
+                if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j))
+                    @. ᶜa_scalar[colidx] =
+                        ᶜh_totʲs.:($$j)[colidx] * draft_area(
+                            Y.c.sgsʲs.:($$j).ρa[colidx],
+                            ᶜρʲs.:($$j)[colidx],
+                        )
+                    vertical_transport!(
+                        Yₜ.c.sgsʲs.:($j).ρae_tot[colidx],
+                        ᶜJ[colidx],
+                        ᶜρʲs.:($j)[colidx],
+                        ᶠu³ʲs.:($j)[colidx],
+                        ᶜa_scalar[colidx],
+                        dt,
+                        edmfx_upwinding,
+                    )
+                end
+
+                for (ᶜρaχʲₜ, ᶜχʲ, χ_name) in
+                    matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j))
+                    χ_name == :e_tot && continue
+                    @. ᶜa_scalar[colidx] =
+                        ᶜχʲ[colidx] * draft_area(
+                            Y.c.sgsʲs.:($$j).ρa[colidx],
+                            ᶜρʲs.:($$j)[colidx],
+                        )
+                    vertical_transport!(
+                        ᶜρaχʲₜ[colidx],
+                        ᶜJ[colidx],
+                        ᶜρʲs.:($j)[colidx],
+                        ᶠu³ʲs.:($j)[colidx],
+                        ᶜa_scalar[colidx],
+                        dt,
+                        edmfx_upwinding,
+                    )
+                end
             end
+        end
 
-            for (ᶜρaχʲₜ, ᶜχʲ, χ_name) in
-                matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j))
-                χ_name == :e_tot && continue
+        if p.atmos.turbconv_model isa AdvectiveEDMFX
+            for j in 1:n
                 @. ᶜa_scalar[colidx] =
-                    ᶜχʲ[colidx] *
                     draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx])
                 vertical_transport!(
-                    ᶜρaχʲₜ[colidx],
+                    Yₜ.c.sgsʲs.:($j).ρa[colidx],
                     ᶜJ[colidx],
                     ᶜρʲs.:($j)[colidx],
                     ᶠu³ʲs.:($j)[colidx],
@@ -182,6 +226,20 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     dt,
                     edmfx_upwinding,
                 )
+
+                vertical_advection!(
+                    Yₜ.c.sgsʲs.:($j).h_tot[colidx],
+                    ᶠu³ʲs.:($j)[colidx],
+                    Y.c.sgsʲs.:($j).h_tot[colidx],
+                    edmfx_upwinding,
+                )
+
+                vertical_advection!(
+                    Yₜ.c.sgsʲs.:($j).q_tot[colidx],
+                    ᶠu³ʲs.:($j)[colidx],
+                    Y.c.sgsʲs.:($j).q_tot[colidx],
+                    edmfx_upwinding,
+                )
             end
         end
 
diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
index 31d332414b7..ffe83940465 100644
--- a/src/prognostic_equations/edmfx_closures.jl
+++ b/src/prognostic_equations/edmfx_closures.jl
@@ -76,10 +76,17 @@ end
 edmfx_nh_pressure_cache(Y, turbconv_model) = (;)
 
 edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
-function edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
+function edmfx_nh_pressure_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::Union{EDMFX, AdvectiveEDMFX},
+)
 
     n = n_mass_flux_subdomains(turbconv_model)
-    (; params, ᶜρʲs, ᶜρ_ref, ᶜp, ᶜp_ref, ᶠgradᵥ_ᶜΦ, ᶠu₃⁰) = p
+    (; params, ᶜρʲs, ᶠgradᵥ_ᶜΦ, ᶠu₃⁰) = p
     FT = eltype(Y)
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
index 44ceaf99b58..a5d9a7317f4 100644
--- a/src/prognostic_equations/edmfx_entr_detr.jl
+++ b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -211,7 +211,7 @@ function detrainment(
         g = CAP.grav(params)
         turbconv_params = CAP.turbconv_params(params)
         ᶜaʲ_max = TCP.max_area(turbconv_params)
-        max_area_limiter = 0.1 * exp(-10 * (ᶜaʲ_max - ᶜaʲ))
+        max_area_limiter = FT(0.1) * exp(-10 * (ᶜaʲ_max - ᶜaʲ))
 
         # pressure scale height (height where pressure drops by 1/e)
         ref_H = ᶜp / (ᶜρ * g)
@@ -254,10 +254,17 @@ function detrainment(
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
     dt::FT,
-    ::ConstantCoefficientDetrainment,
+    ::BOverWDetrainment,
 ) where {FT}
     detr_coeff = CAP.detr_coeff(params)
-    detr = min(detr_coeff * abs(ᶜwʲ), 1 / dt)
+    detr = min(
+        abs(ᶜwʲ) * (
+            detr_coeff +
+            FT(0.12) * abs(min(ᶜbuoyʲ - ᶜbuoy⁰, 0)) /
+            max(eps(FT), (ᶜwʲ - ᶜw⁰) * (ᶜwʲ - ᶜw⁰))
+        ),
+        1 / dt,
+    )
     return detr
 end
 
@@ -283,6 +290,28 @@ function detrainment(
     return detr * FT(2) * hm_limiter(ᶜaʲ)
 end
 
+function detrainment(
+    params,
+    ᶜz::FT,
+    z_sfc::FT,
+    ᶜp::FT,
+    ᶜρ::FT,
+    buoy_flux_surface::FT,
+    ᶜaʲ::FT,
+    ᶜwʲ::FT,
+    ᶜRHʲ::FT,
+    ᶜbuoyʲ::FT,
+    ᶜw⁰::FT,
+    ᶜRH⁰::FT,
+    ᶜbuoy⁰::FT,
+    dt::FT,
+    ::ConstantTimescaleDetrainment,
+) where {FT}
+    detr_tau = CAP.detr_tau(params)
+    detr = min(1 / detr_tau, 1 / dt)
+    return detr
+end
+
 edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
 function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
 
@@ -315,24 +344,36 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
     return nothing
 end
 
-function detrainment(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::ConstantTimescaleDetrainment,
-) where {FT}
-    detr_tau = CAP.detr_tau(params)
-    detr = min(1 / detr_tau, 1 / dt)
-    return detr
+function edmfx_entr_detr_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::AdvectiveEDMFX,
+)
+
+    n = n_mass_flux_subdomains(turbconv_model)
+    (; ᶜentrʲs, ᶜdetrʲs) = p
+    (; ᶜq_tot⁰, ᶜh_tot⁰, ᶠu₃⁰) = p
+
+    for j in 1:n
+
+        @. Yₜ.c.sgsʲs.:($$j).ρa[colidx] +=
+            Y.c.sgsʲs.:($$j).ρa[colidx] *
+            (ᶜentrʲs.:($$j)[colidx] - ᶜdetrʲs.:($$j)[colidx])
+
+        @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] +=
+            ᶜentrʲs.:($$j)[colidx] *
+            (ᶜh_tot⁰[colidx] - Y.c.sgsʲs.:($$j).h_tot[colidx])
+
+        @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] +=
+            ᶜentrʲs.:($$j)[colidx] *
+            (ᶜq_tot⁰[colidx] - Y.c.sgsʲs.:($$j).q_tot[colidx])
+
+        @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] +=
+            ᶠinterp(ᶜentrʲs.:($$j)[colidx]) *
+            (ᶠu₃⁰[colidx] - Y.f.sgsʲs.:($$j).u₃[colidx])
+    end
+    return nothing
 end
diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
index 1034f362448..c1c1d7cc822 100644
--- a/src/prognostic_equations/edmfx_sgs_flux.jl
+++ b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -10,19 +10,18 @@ function edmfx_sgs_mass_flux_tendency!(
     p,
     t,
     colidx,
-    turbconv_model::EDMFX,
+    turbconv_model::Union{EDMFX, AdvectiveEDMFX},
 )
 
     FT = Spaces.undertype(axes(Y.c))
     n = n_mass_flux_subdomains(turbconv_model)
-    (; edmfx_upwinding, sfc_conditions) = p
+    (; edmfx_upwinding) = p
     (; ᶠu³, ᶜh_tot, ᶜspecific) = p
-    (; ᶠu³ʲs, ᶜh_totʲs, ᶜspecificʲs) = p
-    (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰, ᶜspecific⁰) = p
-    (; ᶜK_u, ᶜK_h) = p
+    (; ᶠu³ʲs) = p
+    (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰) = p
+    ᶜq_tot⁰ = turbconv_model isa EDMFX ? p.ᶜspecific⁰.q_tot : p.ᶜq_tot⁰
     (; dt) = p.simulation
     ᶜJ = Fields.local_geometry_field(Y.c).J
-    ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_sgs_mass_flux
         # energy
@@ -30,7 +29,10 @@ function edmfx_sgs_mass_flux_tendency!(
         ᶜh_tot_diff_colidx = ᶜq_tot_diff_colidx = p.ᶜtemp_scalar[colidx]
         for j in 1:n
             @. ᶠu³_diff_colidx = ᶠu³ʲs.:($$j)[colidx] - ᶠu³[colidx]
-            @. ᶜh_tot_diff_colidx = ᶜh_totʲs.:($$j)[colidx] - ᶜh_tot[colidx]
+            ᶜh_totʲ_colidx =
+                turbconv_model isa EDMFX ? p.ᶜh_totʲs.:($j)[colidx] :
+                Y.c.sgsʲs.:($j).h_tot[colidx]
+            @. ᶜh_tot_diff_colidx = ᶜh_totʲ_colidx - ᶜh_tot[colidx]
             vertical_transport!(
                 Yₜ.c.ρe_tot[colidx],
                 ᶜJ[colidx],
@@ -57,8 +59,11 @@ function edmfx_sgs_mass_flux_tendency!(
             # specific humidity
             for j in 1:n
                 @. ᶠu³_diff_colidx = ᶠu³ʲs.:($$j)[colidx] - ᶠu³[colidx]
-                @. ᶜq_tot_diff_colidx =
-                    ᶜspecificʲs.:($$j).q_tot[colidx] - ᶜspecific.q_tot[colidx]
+                ᶜq_totʲ_colidx =
+                    turbconv_model isa EDMFX ?
+                    p.ᶜspecificʲs.:($j).q_tot[colidx] :
+                    Y.c.sgsʲs.:($j).q_tot[colidx]
+                @. ᶜq_tot_diff_colidx = ᶜq_totʲ_colidx - ᶜspecific.q_tot[colidx]
                 vertical_transport!(
                     Yₜ.c.ρq_tot[colidx],
                     ᶜJ[colidx],
@@ -70,8 +75,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 )
             end
             @. ᶠu³_diff_colidx = ᶠu³⁰[colidx] - ᶠu³[colidx]
-            @. ᶜq_tot_diff_colidx =
-                ᶜspecific⁰.q_tot[colidx] - ᶜspecific.q_tot[colidx]
+            @. ᶜq_tot_diff_colidx = ᶜq_tot⁰[colidx] - ᶜspecific.q_tot[colidx]
             vertical_transport!(
                 Yₜ.c.ρq_tot[colidx],
                 ᶜJ[colidx],
@@ -159,18 +163,14 @@ function edmfx_sgs_diffusive_flux_tendency!(
     p,
     t,
     colidx,
-    turbconv_model::EDMFX,
+    turbconv_model::Union{EDMFX, AdvectiveEDMFX},
 )
 
     FT = Spaces.undertype(axes(Y.c))
-    n = n_mass_flux_subdomains(turbconv_model)
-    (; edmfx_upwinding, sfc_conditions) = p
-    (; ᶠu³, ᶜh_tot, ᶜspecific) = p
-    (; ᶠu³ʲs, ᶜh_totʲs, ᶜspecificʲs) = p
-    (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰, ᶜspecific⁰) = p
+    (; sfc_conditions) = p
+    (; ᶜρa⁰, ᶠu³⁰, ᶜu⁰, ᶜh_tot⁰) = p
+    ᶜq_tot⁰ = turbconv_model isa EDMFX ? p.ᶜspecific⁰.q_tot : p.ᶜq_tot⁰
     (; ᶜK_u, ᶜK_h) = p
-    (; dt) = p.simulation
-    ᶜJ = Fields.local_geometry_field(Y.c).J
     ᶠgradᵥ = Operators.GradientC2F()
 
     if p.atmos.edmfx_sgs_diffusive_flux
@@ -193,9 +193,8 @@ function edmfx_sgs_diffusive_flux_tendency!(
                     sfc_conditions.ρ_flux_q_tot[colidx],
                 ),
             )
-            @. Yₜ.c.ρq_tot[colidx] -= ᶜdivᵥ_ρq_tot(
-                -(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜspecific⁰.q_tot[colidx])),
-            )
+            @. Yₜ.c.ρq_tot[colidx] -=
+                ᶜdivᵥ_ρq_tot(-(ᶠρaK_h[colidx] * ᶠgradᵥ(ᶜq_tot⁰[colidx])))
         end
 
         # momentum
diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
index c7a886184d6..5ff7aca7dab 100644
--- a/src/prognostic_equations/hyperdiffusion.jl
+++ b/src/prognostic_equations/hyperdiffusion.jl
@@ -21,7 +21,10 @@ function hyperdiffusion_cache(Y, atmos, do_dss)
 
     # Sub-grid scale quantities
     ᶜ∇²uʲs =
-        atmos.turbconv_model isa EDMFX ? similar(Y.c, NTuple{n, C123{FT}}) : (;)
+        (
+            atmos.turbconv_model isa EDMFX ||
+            atmos.turbconv_model isa AdvectiveEDMFX
+        ) ? similar(Y.c, NTuple{n, C123{FT}}) : (;)
     sgs_quantities =
         atmos.turbconv_model isa EDMFX ?
         (;
@@ -33,6 +36,14 @@ function hyperdiffusion_cache(Y, atmos, do_dss)
                 specific_sgsʲs.(Y.c, atmos.turbconv_model)
             ),
         ) :
+        atmos.turbconv_model isa AdvectiveEDMFX ?
+        (;
+            ᶜ∇²tke⁰ = similar(Y.c, FT),
+            ᶜ∇²uₕʲs = similar(Y.c, NTuple{n, C12{FT}}),
+            ᶜ∇²uᵥʲs = similar(Y.c, NTuple{n, C3{FT}}),
+            ᶜ∇²h_totʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜ∇²q_totʲs = similar(Y.c, NTuple{n, FT}),
+        ) :
         atmos.turbconv_model isa DiagnosticEDMFX ?
         (; ᶜ∇²tke⁰ = similar(Y.c, FT)) : (;)
     quantities = (; gs_quantities..., sgs_quantities...)
@@ -71,6 +82,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
             ᶜ∇²specific_energyʲs,
         ) = p
     end
+    if turbconv_model isa AdvectiveEDMFX
+        (; ᶜρa⁰, ᶜρʲs, ᶜ∇²tke⁰, ᶜtke⁰, ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²h_totʲs) = p
+    end
     if turbconv_model isa DiagnosticEDMFX
         (; ᶜtke⁰, ᶜ∇²tke⁰) = p
     end
@@ -89,7 +103,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     if turbconv_model isa EDMFX && diffuse_tke
         @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜspecific⁰.tke))
     end
-    if turbconv_model isa DiagnosticEDMFX && diffuse_tke
+    if (
+        turbconv_model isa DiagnosticEDMFX || turbconv_model isa AdvectiveEDMFX
+    ) && diffuse_tke
         @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜtke⁰))
     end
 
@@ -110,6 +126,15 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
         end
     end
 
+    if turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            @. ᶜ∇²uʲs.:($$j) =
+                C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j)))) -
+                C123(wcurlₕ(C123(curlₕ(p.ᶜuʲs.:($$j)))))
+            @. ᶜ∇²h_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).h_tot))
+        end
+    end
+
     if do_dss
         NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin
             for dss_op! in (
@@ -122,10 +147,11 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
                 # operations do not accept Covariant123Vector types
                 dss_op!(ᶜ∇²u, buffer.ᶜ∇²u)
                 dss_op!(ᶜ∇²specific_energy, buffer.ᶜ∇²specific_energy)
-                if turbconv_model isa EDMFX && diffuse_tke
+                if diffuse_tke
                     dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰)
                 end
-                if turbconv_model isa EDMFX
+                if (turbconv_model isa EDMFX) ||
+                   (turbconv_model isa AdvectiveEDMFX)
                     # Need to split the DSS computation here, because our DSS 
                     # operations do not accept Covariant123Vector types     
                     for j in 1:n
@@ -138,11 +164,16 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
                         @. ᶜ∇²uʲs.:($$j) =
                             C123(ᶜ∇²uₕʲs.:($$j)) + C123(ᶜ∇²uᵥʲs.:($$j))
                     end
+                end
+                if turbconv_model isa EDMFX
                     dss_op!(ᶜ∇²specific_energyʲs, buffer.ᶜ∇²specific_energyʲs)
                 end
-                if turbconv_model isa DiagnosticEDMFX && diffuse_tke
-                    dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰)
+                if turbconv_model isa AdvectiveEDMFX
+                    dss_op!(ᶜ∇²h_totʲs, buffer.ᶜ∇²h_totʲs)
                 end
+                # if turbconv_model isa DiagnosticEDMFX && diffuse_tke
+                #     dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰)
+                # end
             end
         end
     end
@@ -158,7 +189,8 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
     @. ᶜρ_energyₜ -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy))
 
     # Sub-grid scale hyperdiffusion continued
-    if turbconv_model isa EDMFX && diffuse_tke
+    if (turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX) &&
+       diffuse_tke
         @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(ᶜρa⁰ * gradₕ(ᶜ∇²tke⁰))
     end
     if turbconv_model isa EDMFX
@@ -177,6 +209,19 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
                 wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²specific_energyʲs.:($$j)))
         end
     end
+    if turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            if point_type <: Geometry.Abstract3DPoint
+                # only need curl-curl part
+                @. ᶜ∇²uᵥʲs.:($$j) = C3(wcurlₕ(C123(curlₕ(ᶜ∇²uʲs.:($$j)))))
+                @. Yₜ.f.sgsʲs.:($$j).u₃ +=
+                    κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, ᶜ∇²uᵥʲs.:($$j))
+            end
+            # Note: It is more correct to have ρa inside and outside the divergence
+            @. Yₜ.c.sgsʲs.:($$j).h_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²h_totʲs.:($$j)))
+        end
+    end
+
     if turbconv_model isa DiagnosticEDMFX && diffuse_tke
         @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰))
     end
@@ -193,6 +238,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     if turbconv_model isa EDMFX
         (; ᶜspecificʲs, ᶜ∇²specific_tracersʲs, ᶜp) = p
     end
+    if turbconv_model isa AdvectiveEDMFX
+        (; ᶜ∇²q_totʲs) = p
+    end
     if do_dss
         buffer = p.hyperdiffusion_ghost_buffer
     end
@@ -209,6 +257,13 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
         end
     end
 
+    if turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            # Note: It is more correct to have ρa inside and outside the divergence
+            @. ᶜ∇²q_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_tot))
+        end
+    end
+
     if do_dss
         NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin
             for dss_op! in (
@@ -220,6 +275,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
                 if turbconv_model isa EDMFX
                     dss_op!(ᶜ∇²specific_tracersʲs, buffer.ᶜ∇²specific_tracersʲs)
                 end
+                if turbconv_model isa AdvectiveEDMFX
+                    dss_op!(ᶜ∇²q_totʲs, buffer.ᶜ∇²q_totʲs)
+                end
             end
         end
     end
@@ -245,5 +303,12 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
             end
         end
     end
+    if turbconv_model isa AdvectiveEDMFX
+        for j in 1:n
+            @. Yₜ.c.sgsʲs.:($$j).ρa -=
+                κ₄ * wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²q_totʲs.:($$j)))
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j)))
+        end
+    end
     return nothing
 end
diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl
index c7c97a68f56..ca89f3582b5 100644
--- a/src/prognostic_equations/implicit/implicit_tendency.jl
+++ b/src/prognostic_equations/implicit/implicit_tendency.jl
@@ -52,6 +52,13 @@ vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak})
         ),
     ))
 
+vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:none}) =
+    @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠu³ * ᶠinterp(ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)
+vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:first_order}) =
+    @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠupwind1(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)
+vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:third_order}) =
+    @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠupwind3(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)
+
 function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     (; energy_upwinding, tracer_upwinding, density_upwinding, edmfx_upwinding) =
         p
@@ -107,7 +114,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     if rayleigh_sponge isa RayleighSponge
         (; ᶠβ_rayleigh_w) = p
         @. Yₜ.f.u₃[colidx] -= ᶠβ_rayleigh_w[colidx] * Y.f.u₃[colidx]
-        if turbconv_model isa EDMFX
+        if turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX
             for j in 1:n
                 @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
                     ᶠβ_rayleigh_w[colidx] * Y.f.sgsʲs.:($$j).u₃[colidx]
diff --git a/src/prognostic_equations/pressure_work.jl b/src/prognostic_equations/pressure_work.jl
index a190c2cc5f9..6df349c1e0f 100644
--- a/src/prognostic_equations/pressure_work.jl
+++ b/src/prognostic_equations/pressure_work.jl
@@ -14,4 +14,7 @@ function pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::EDMFX)
     return nothing
 end
 
+# We ignore the small pressure term in advective EDMF for now
+pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::AdvectiveEDMFX) = nothing
+
 pressure_work_tendency!(Yₜ, Y, p, t, colidx, ::Any) = nothing
diff --git a/src/prognostic_equations/zero_velocity.jl b/src/prognostic_equations/zero_velocity.jl
index f62a79b55cc..cb6787ddeb1 100644
--- a/src/prognostic_equations/zero_velocity.jl
+++ b/src/prognostic_equations/zero_velocity.jl
@@ -10,7 +10,8 @@ function zero_velocity_tendency!(Yₜ, Y, p, t, colidx)
 
         @. Yₜ.c.uₕ[colidx] = C12(FT(0), FT(0))
         @. Yₜ.f.u₃[colidx] = Geometry.Covariant3Vector(FT(0))
-        if p.atmos.turbconv_model isa EDMFX
+        if p.atmos.turbconv_model isa EDMFX ||
+           p.atmos.turbconv_model isa AdvectiveEDMFX
             for j in 1:n
                 @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] =
                     Geometry.Covariant3Vector(FT(0))
diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
index ca0600dfb28..e1c0d7f35e6 100644
--- a/src/solver/model_getters.jl
+++ b/src/solver/model_getters.jl
@@ -341,7 +341,8 @@ function get_turbconv_model(
     turbconv_params,
 )
     turbconv = parsed_args["turbconv"]
-    @assert turbconv in (nothing, "edmf", "edmfx", "diagnostic_edmfx")
+    @assert turbconv in
+            (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx")
 
     return if turbconv == "edmf"
         TC.EDMFModel(
@@ -355,6 +356,10 @@ function get_turbconv_model(
         N = turbconv_params.updraft_number
         TKE = parsed_args["prognostic_tke"]
         EDMFX{N, TKE}(turbconv_params.min_area)
+    elseif turbconv == "advective_edmfx"
+        N = turbconv_params.updraft_number
+        TKE = parsed_args["prognostic_tke"]
+        AdvectiveEDMFX{N, TKE}(turbconv_params.min_area)
     elseif turbconv == "diagnostic_edmfx"
         N = turbconv_params.updraft_number
         TKE = parsed_args["prognostic_tke"]
@@ -387,14 +392,14 @@ function get_detrainment_model(parsed_args)
         NoDetrainment()
     elseif detr_model == "PiGroups"
         PiGroupsDetrainment()
-    elseif detr_model == "ConstantCoefficient"
-        ConstantCoefficientDetrainment()
+    elseif detr_model == "BOverW"
+        BOverWDetrainment()
     elseif detr_model == "ConstantCoefficientHarmonics"
         ConstantCoefficientHarmonicsDetrainment()
     elseif detr_model == "ConstantTimescale"
         ConstantTimescaleDetrainment()
     else
-        error("Invalid entr_model $(entr_model)")
+        error("Invalid detr_model $(detr_model)")
     end
 end
 
diff --git a/src/solver/types.jl b/src/solver/types.jl
index 9673374188e..0d4ee597941 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -157,6 +157,12 @@ struct EDMFX{N, TKE, FT} <: AbstractEDMF
 end
 EDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = EDMFX{N, TKE, FT}(a_half)
 
+struct AdvectiveEDMFX{N, TKE, FT} <: AbstractEDMF
+    a_half::FT # WARNING: this should never be used outside of divide_by_ρa
+end
+AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} =
+    AdvectiveEDMFX{N, TKE, FT}(a_half)
+
 struct DiagnosticEDMFX{N, TKE, FT} <: AbstractEDMF
     a_int::FT # area fraction of the first interior cell above the surface
     a_half::FT # WARNING: this should never be used outside of divide_by_ρa
@@ -165,13 +171,16 @@ DiagnosticEDMFX{N, TKE}(a_int::FT, a_half::FT) where {N, TKE, FT} =
     DiagnosticEDMFX{N, TKE, FT}(a_int, a_half)
 
 n_mass_flux_subdomains(::EDMFX{N}) where {N} = N
+n_mass_flux_subdomains(::AdvectiveEDMFX{N}) where {N} = N
 n_mass_flux_subdomains(::DiagnosticEDMFX{N}) where {N} = N
 n_mass_flux_subdomains(::Any) = 0
 
 n_prognostic_mass_flux_subdomains(::EDMFX{N}) where {N} = N
+n_prognostic_mass_flux_subdomains(::AdvectiveEDMFX{N}) where {N} = N
 n_prognostic_mass_flux_subdomains(::Any) = 0
 
 use_prognostic_tke(::EDMFX{N, TKE}) where {N, TKE} = TKE
+use_prognostic_tke(::AdvectiveEDMFX{N, TKE}) where {N, TKE} = TKE
 use_prognostic_tke(::DiagnosticEDMFX{N, TKE}) where {N, TKE} = TKE
 use_prognostic_tke(::Any) = false
 
@@ -186,7 +195,7 @@ abstract type AbstractDetrainmentModel end
 
 struct NoDetrainment <: AbstractDetrainmentModel end
 struct PiGroupsDetrainment <: AbstractDetrainmentModel end
-struct ConstantCoefficientDetrainment <: AbstractDetrainmentModel end
+struct BOverWDetrainment <: AbstractDetrainmentModel end
 struct ConstantCoefficientHarmonicsDetrainment <: AbstractDetrainmentModel end
 struct ConstantTimescaleDetrainment <: AbstractDetrainmentModel end
 
@@ -234,6 +243,7 @@ Base.broadcastable(x::AbstractEnergyFormulation) = tuple(x)
 Base.broadcastable(x::AbstractPrecipitationModel) = tuple(x)
 Base.broadcastable(x::AbstractForcing) = tuple(x)
 Base.broadcastable(x::EDMFX) = tuple(x)
+Base.broadcastable(x::AdvectiveEDMFX) = tuple(x)
 Base.broadcastable(x::DiagnosticEDMFX) = tuple(x)
 Base.broadcastable(x::AbstractEntrainmentModel) = tuple(x)
 Base.broadcastable(x::AbstractDetrainmentModel) = tuple(x)
diff --git a/toml/advective_edmfx_box.toml b/toml/advective_edmfx_box.toml
new file mode 100644
index 00000000000..3c8a8e76d5d
--- /dev/null
+++ b/toml/advective_edmfx_box.toml
@@ -0,0 +1,24 @@
+[C_E]
+alias = "C_E"
+value = 0.044
+type = "float"
+
+[EDMF_surface_area]
+alias = "surface_area"
+value = 0.1
+type = "float"
+
+[EDMF_min_area]
+alias = "min_area"
+value = 1.0e-3
+type = "float"
+
+[entr_coeff]
+alias = "entr_coeff"
+value = 0.1
+type = "float"
+
+[entr_tau]
+alias = "entr_tau"
+value = 1800
+type = "float"