diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
index f721db6ab61..ca884961772 100644
--- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
+++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl
@@ -165,6 +165,8 @@ function vars_state_auxiliary(m::HBModel, T)
         w::T     # ∫(-∇⋅u)
         pkin::T  # ∫(-αᵀθ)
         wz0::T   # w at z=0
+        uᵈ::SVector{2, T}    # velocity deviation from vertical mean
+        ΔGᵘ::SVector{2, T}   # vertically averaged tendency
     end
 end
 
diff --git a/src/Ocean/ShallowWater/ShallowWaterModel.jl b/src/Ocean/ShallowWater/ShallowWaterModel.jl
index 7616c5a34aa..462e81bb790 100644
--- a/src/Ocean/ShallowWater/ShallowWaterModel.jl
+++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl
@@ -56,8 +56,8 @@ SWModel = ShallowWaterModel
 
 function vars_state_conservative(m::SWModel, T)
     @vars begin
-        U::SVector{3, T}
         η::T
+        U::SVector{3, T}
     end
 end
 
@@ -65,6 +65,12 @@ function vars_state_auxiliary(m::SWModel, T)
     @vars begin
         f::SVector{3, T}
         τ::SVector{3, T}  # value includes τₒ, g, and ρ
+        Gᵁ::SVector{3, T} # integral of baroclinic tendency
+        ηᶜ::T             # cummulate η value over fast steps
+        Uᶜ::SVector{3, T} # cummulate U value over fast steps
+        ηˢ::T             # starting η field value
+        Uˢ::SVector{3, T} # starting U field value
+        Δu::SVector{3, T} # reconciliation Δu = 1/H * (Ū - ∫u)
     end
 end
 
diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
index cbf05ee265e..116d85c423a 100644
--- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl
+++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl
@@ -9,28 +9,134 @@ import ..DGMethods:
     cummulate_fast_solution!,
     reconcile_from_fast_to_slow!
 
-@inline initialize_states!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
-    _...,
-) = nothing
-
-@inline tendency_from_slow_to_fast!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
-    _...,
-) = nothing
-
-@inline cummulate_fast_solution!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
-    _...,
-) = nothing
-
-@inline reconcile_from_fast_to_slow!(
-    slow::HydrostaticBoussinesqModel,
-    fast::ShallowWaterModel,
-    _...,
-) = nothing
+@inline function initialize_states!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+)
+    model_bt.state_auxiliary.ηᶜ .= -0
+    model_bt.state_auxiliary.Uᶜ .= (@SVector [-0, -0, -0])'
+
+    state_bt.η .= model_bt.state_auxiliary.ηˢ
+    state_bt.U .= model_bt.state_auxiliary.Uˢ
+
+    model_bc.state_auxiliary.ΔGᵘ .= -0
+
+    return nothing
+end
+
+@inline function tendency_from_slow_to_fast!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+    forcing_tendency,
+)
+    # integrate the tendency
+    tendency_model = model_bc.modeldata.tendency_model
+    update_auxiliary_state!(
+        tendency_model,
+        tendency_model.balance_law,
+        forcing_tendency,
+        0,
+    )
+
+    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
+
+    ## get top value (=integral over full depth) of ∫du
+    ∫du = tendency_model.state_auxiliary.∫du
+    boxy_∫du = reshape(∫du, Nq^2, Nqk, 2, nelemv, nelemh)
+    flat_∫du = @view boxy_∫du[:, end, :, end, :]
+
+    ## copy into Gᵁ of dgFast
+    model_bt.state_auxiliary.Gᵁ .= reshape(flat_∫du, Nq^2, 2, nelemh)
+
+    ## scale by -1/H and copy back to ΔGu
+    # note: since tendency_dg.auxstate.∫du is not used after this, could be
+    # re-used to store a 3-D copy of "-Gu"
+    ΔGᵘ = model_bc.state_auxiliary.ΔGᵘ
+    boxy_ΔGᵘ = reshape(ΔGᵘ, Nq^2, Nqk, 2, nelemv, nelemh)
+    boxy_ΔGᵘ .= -reshape(flat_∫du, Nq^2, 1, 2, 1, nelemh) / baroclinic.problem.H
+
+    return nothing
+end
+
+@inline function cummulate_fast_solution!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bt,
+    state_bt,
+    fast_time,
+    fast_dt,
+    substep,
+)
+    ### for now no averaging, just save the final time step
+    model_bt.state_auxiliary.ηᶜ .= state_bt.η
+    model_bt.state_auxiliary.Uᶜ .= state_bt.U
+
+    ### technically unnecessary right now, but for future changes
+    model_bt.state_auxiliary.ηˢ .= state_bt.η
+    model_bt.state_auxiliary.Uˢ .= state_bt.U
+
+    return nothing
+end
+
+@inline function reconcile_from_fast_to_slow!(
+    baroclinic::HydrostaticBoussinesqModel,
+    barotropic::ShallowWaterModel,
+    model_bc,
+    model_bt,
+    state_bc,
+    state_bt,
+)
+    Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(model_bc)
+
+    # need to calculate int_u using integral kernels
+    # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow)
+
+    # Compute: \int_{-H}^{0} u_slow)
+    ### need to make sure this is stored into aux.∫u
+
+    # integrate vertically horizontal velocity
+    flow_model = model_bc.modeldata.flow_model
+    update_auxiliary_state!(flow_model, flow_model.balance_law, state_bc, 0)
+
+    ## get top value (=integral over full depth)
+    ∫u = flow_model.state_auxiliary.∫u
+    boxy_∫u = reshape(∫u, Nq^2, Nqk, 2, nelemv, nelemh)
+    flat_∫u = @view boxy_∫u[:, end, :, end, :]
+
+    ## get time weighted averaged out of cumulative arrays
+    ## no weights right now
+    state_bt.state_auxiliary.Uᶜ .*= 1 # / fast_time_rec[2]
+    state_bt.state_auxiliary.ηᶜ .*= 1 # / fast_time_rec[2]
+
+    ### substract ∫u from U and divide by H
+
+    ### Δu is a place holder for 1/H * (Ū - ∫u)
+    Δu = model_bt.state_auxiliary.Δu
+    Δu .= 1 / baroclinic.problem.H * (model_bt.state_auxiliary.Uᶜ - flat_∫u)
+
+    ### copy the 2D contribution down the 3D solution
+    ### need to reshape these things for the broadcast
+    boxy_u = reshape(state_bc.u, Nq^2, Nqk, 2, nelemv, nelemh)
+    boxy_Δu = reshape(Δu, Nq^2, 1, 3, 1, nelemh)
+    ### this works, we tested it
+    # currently SWModel has a 3-vector for U
+    # need to shrink to 2-vector at some point
+    boxy_u .+= boxy_Δu[:, :, 1:2, :, :]
+
+    ### copy 2D eta over to 3D model
+    boxy_η_3D = reshape(state_bc.η, Nq^2, Nqk, nelemv, nelemh)
+    boxy_η_2D = reshape(model_bt.state_auxiliary.ηᶜ, Nq^2, 1, 1, nelemh)
+    boxy_η_3D .= boxy_η_2D
+
+    return nothing
+end
 
 end