diff --git a/docs/src/APIs/Numerics/DGMethods/DGMethods.md b/docs/src/APIs/Numerics/DGMethods/DGMethods.md
index 74940efd62e..1ee0da3939f 100644
--- a/docs/src/APIs/Numerics/DGMethods/DGMethods.md
+++ b/docs/src/APIs/Numerics/DGMethods/DGMethods.md
@@ -5,6 +5,7 @@ CurrentModule = ClimateMachine.DGMethods
 ```
 
 ```@docs
+DGModel
 remainder_DGModel
 ```
 
diff --git a/src/Numerics/DGMethods/DGModel.jl b/src/Numerics/DGMethods/DGModel.jl
index ea276cfb7b0..b45bc70f92a 100644
--- a/src/Numerics/DGMethods/DGModel.jl
+++ b/src/Numerics/DGMethods/DGModel.jl
@@ -46,13 +46,31 @@ end
 # Include the remainder model for composing DG models and balance laws
 include("remainder.jl")
 
+"""
+    (dg::DGModel)(tendency, state_conservative, nothing, t, α, β)
+
+Computes the tendency terms compatible with `IncrementODEProblem`
+
+    tendency .= α .* dQdt(state_conservative, p, t) .+ β .* tendency
+
+The 4-argument form will just compute
+
+    tendency .= dQdt(state_conservative, p, t) 
+
+"""
 function (dg::DGModel)(
     tendency,
     state_conservative,
-    ::Nothing,
+    param,
     t;
     increment = false,
 )
+    # TODO deprecate increment argument
+    dg(tendency, state_conservative, param, t, true, increment)
+end
+
+function (dg::DGModel)(tendency, state_conservative, _, t, α, β)
+
 
     balance_law = dg.balance_law
     device = array_device(state_conservative)
@@ -87,8 +105,10 @@ function (dg::DGModel)(
     ndrange_interior_surface = Nfp * length(grid.interiorelems)
     ndrange_exterior_surface = Nfp * length(grid.exteriorelems)
 
-    if !increment && num_state_conservative < num_state_tendency
-        tendency .= -zero(FT)
+    if num_state_conservative < num_state_tendency && β != 1
+        # if we don't operate on the full state, then we need to scale here instead of volume_tendency!
+        tendency .*= β
+        β = β != 0 # if β==0 then we can avoid the memory load in volume_tendency!  
     end
 
     communicate =
@@ -412,7 +432,8 @@ function (dg::DGModel)(
         grid.ω,
         grid.D,
         topology.realelems,
-        increment;
+        α,
+        β;
         ndrange = (nrealelem * Nq, Nq),
         dependencies = (comp_stream,),
     )
@@ -435,7 +456,8 @@ function (dg::DGModel)(
         grid.vmap⁻,
         grid.vmap⁺,
         grid.elemtobndy,
-        grid.interiorelems;
+        grid.interiorelems,
+        α;
         ndrange = ndrange_interior_surface,
         dependencies = (comp_stream,),
     )
@@ -506,7 +528,8 @@ function (dg::DGModel)(
         grid.vmap⁻,
         grid.vmap⁺,
         grid.elemtobndy,
-        grid.exteriorelems;
+        grid.exteriorelems,
+        α;
         ndrange = ndrange_exterior_surface,
         dependencies = (
             comp_stream,
diff --git a/src/Numerics/DGMethods/DGModel_kernels.jl b/src/Numerics/DGMethods/DGModel_kernels.jl
index 83b2e44c183..d7e137214f7 100644
--- a/src/Numerics/DGMethods/DGModel_kernels.jl
+++ b/src/Numerics/DGMethods/DGModel_kernels.jl
@@ -47,7 +47,8 @@ Computational kernel: Evaluate the volume integrals on right-hand side of a
     ω,
     D,
     elems,
-    increment,
+    α,
+    β,
 ) where {dim, polyorder}
     @uniform begin
         N = polyorder
@@ -88,8 +89,7 @@ Computational kernel: Evaluate the volume integrals on right-hand side of a
         @unroll for k in 1:Nqk
             ijk = i + Nq * ((j - 1) + Nq * (k - 1))
             @unroll for s in 1:num_state_conservative
-                local_tendency[k, s] =
-                    increment ? tendency[ijk, s, e] : zero(FT)
+                local_tendency[k, s] = zero(FT)
             end
             local_MI[k] = vgeo[ijk, _MI, e]
         end
@@ -302,7 +302,12 @@ Computational kernel: Evaluate the volume integrals on right-hand side of a
         @unroll for k in 1:Nqk
             ijk = i + Nq * ((j - 1) + Nq * (k - 1))
             @unroll for s in 1:num_state_conservative
-                tendency[ijk, s, e] = local_tendency[k, s]
+                if β != 0
+                    T = α * local_tendency[k, s] + β * tendency[ijk, s, e]
+                else
+                    T = α * local_tendency[k, s]
+                end
+                tendency[ijk, s, e] = T
             end
         end
     end
@@ -323,7 +328,8 @@ end
     ω,
     D,
     elems,
-    increment,
+    α,
+    β,
 ) where {dim, polyorder}
     @uniform begin
         N = polyorder
@@ -372,8 +378,7 @@ end
         @unroll for k in 1:Nqk
             ijk = i + Nq * ((j - 1) + Nq * (k - 1))
             @unroll for s in 1:num_state_conservative
-                local_tendency[k, s] =
-                    increment ? tendency[ijk, s, e] : zero(FT)
+                local_tendency[k, s] = zero(FT)
             end
             local_MI[k] = vgeo[ijk, _MI, e]
         end
@@ -509,7 +514,12 @@ end
         @unroll for k in 1:Nqk
             ijk = i + Nq * ((j - 1) + Nq * (k - 1))
             @unroll for s in 1:num_state_conservative
-                tendency[ijk, s, e] = local_tendency[k, s]
+                if β != 0
+                    T = α * local_tendency[k, s] + β * tendency[ijk, s, e]
+                else
+                    T = α * local_tendency[k, s]
+                end
+                tendency[ijk, s, e] = T
             end
         end
     end
@@ -545,6 +555,7 @@ Computational kernel: Evaluate the surface integrals on right-hand side of a
     vmap⁺,
     elemtobndy,
     elems,
+    α,
 ) where {dim, polyorder}
     @uniform begin
         N = polyorder
@@ -818,7 +829,7 @@ Computational kernel: Evaluate the surface integrals on right-hand side of a
         #Update RHS
         @unroll for s in 1:num_state_conservative
             # FIXME: Should we pretch these?
-            tendency[vid⁻, s, e⁻] -= vMI * sM * local_flux[s]
+            tendency[vid⁻, s, e⁻] -= α * vMI * sM * local_flux[s]
         end
         # Need to wait after even faces to avoid race conditions
         @synchronize(f % 2 == 0)