diff --git a/src/Biogeochemistry.jl b/src/Biogeochemistry.jl
index 21f7e27471..5b0172d7d7 100644
--- a/src/Biogeochemistry.jl
+++ b/src/Biogeochemistry.jl
@@ -59,10 +59,12 @@ is a subtype of `AbstractBioeochemistry`:
 abstract type AbstractBiogeochemistry end
 
 # Returns the forcing for discrete form models
-@inline biogeochemical_transition(i, j, k, grid, bgc, val_tracer_name, clock, fields) =
+@inline biogeochemical_transition(i, j, k, grid, bgc, val_tracer_name, clock, fields, val_compute_bgc) =
     bgc(i, j, k, grid, val_tracer_name, clock, fields)
 
-@inline biogeochemical_transition(i, j, k, grid, ::Nothing, val_tracer_name, clock, fields) = zero(grid)
+@inline biogeochemical_transition(i, j, k, grid, ::Nothing, val_tracer_name, clock, fields, val_compute_bgc) = zero(grid)
+
+@inline biogeochemical_transition(i, j, k, grid, bgc, val_tracer_name, clock, fields, ::Val{false}) = zero(grid)
 
 # Required for when a model is defined but not for all tracers
 @inline (bgc::AbstractBiogeochemistry)(i, j, k, grid, val_tracer_name, clock, fields) = zero(grid)
@@ -106,7 +108,7 @@ abstract type AbstractContinuousFormBiogeochemistry <: AbstractBiogeochemistry e
 
 """Return the biogeochemical forcing for `val_tracer_name` for continuous form when model is called."""
 @inline function biogeochemical_transition(i, j, k, grid, bgc::AbstractContinuousFormBiogeochemistry,
-                                           val_tracer_name, clock, fields)
+                                           val_tracer_name, clock, fields, val_compute_bgc)
 
     names_to_extract = tuple(required_biogeochemical_tracers(bgc)...,
                              required_biogeochemical_auxiliary_fields(bgc)...)
@@ -120,6 +122,8 @@ abstract type AbstractContinuousFormBiogeochemistry <: AbstractBiogeochemistry e
     return bgc(val_tracer_name, x, y, z, clock.time, fields_ijk...)
 end
 
+@inline biogeochemical_transition(i, j, k, grid, bgc::AbstractContinuousFormBiogeochemistry, val_tracer_name, clock, fields, ::Val{false}) = zero(grid)
+
 @inline (bgc::AbstractContinuousFormBiogeochemistry)(val_tracer_name, x, y, z, t, fields...) = zero(t)
 
 tracernames(tracers) = keys(tracers)
@@ -171,4 +175,70 @@ const AbstractBGCOrNothing = Union{Nothing, AbstractBiogeochemistry}
 required_biogeochemical_tracers(::AbstractBGCOrNothing) = ()
 required_biogeochemical_auxiliary_fields(::AbstractBGCOrNothing) = ()
 
+#####
+##### BGC only stepping
+#####
+
+using Oceananigans: interior_tendency_kernel_parameters
+using Oceananigans.Utils: launch!
+
+using KernelAbstractions: @kernel, @index
+
+import Oceananigans: compute_biogeochemical_tendencies!
+
+function compute_biogeochemical_tendencies!(model, tendencies; active_cells_map = nothing)
+
+    arch                 = model.architecture
+    grid                 = model.grid
+    biogeochemistry      = model.biogeochemistry
+    velocities           = model.velocities
+    tracers              = model.tracers
+    auxiliary_fields     = model.auxiliary_fields
+    clock                = model.clock
+
+    kernel_parameters = tuple(interior_tendency_kernel_parameters(grid))
+
+    tracer_kernel_args   = (biogeochemistry, velocities, tracers, auxiliary_fields)
+
+    for tracer_index in 1:length(tracers)
+        @inbounds c_tendency = tendencies[tracer_index + 3]
+        @inbounds tracer_name = keys(tracers)[tracer_index]
+
+        args = tuple(Val(tracer_name), tracer_kernel_args..., clock)
+
+        for parameters in kernel_parameters
+            launch!(arch, grid, parameters, compute_bgc_G!, 
+                    c_tendency, grid, active_cells_map, args;
+                    active_cells_map)
+        end
+    end
+
+    return nothing
+end
+
+""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
+@kernel function compute_bgc_G!(Gc, grid, ::Nothing, args)
+    i, j, k = @index(Global, NTuple)
+    @inbounds Gc[i, j, k] = biogeochemical_tendency(i, j, k, grid, args...)
+end
+
+@kernel function compute_bgc_G!(Gc, grid, interior_map, args) 
+    idx = @index(Global, Linear)
+    i, j, k = active_linear_index_to_tuple(idx, interior_map)
+    @inbounds Gc[i, j, k] = biogeochemical_tendency(i, j, k, grid, args...)
+end
+
+@inline function biogeochemical_tendency(i, j, k, grid,
+                                         val_tracer_name,
+                                         biogeochemistry,
+                                         velocities,
+                                         tracers,
+                                         auxiliary_fields,
+                                         clock)
+
+    model_fields = merge(velocities, tracers, auxiliary_fields)
+
+    return biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields, Val(true))
+end
+
 end # module
diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl
index db53b860d6..8966f63823 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl
@@ -10,6 +10,8 @@ using Oceananigans.Fields: immersed_boundary_condition
 using Oceananigans.Biogeochemistry: update_tendencies!
 using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: FlavorOfCATKE, FlavorOfTD
 
+using Oceananigans.TimeSteppers: compute_bgc_with_physics
+
 using Oceananigans.ImmersedBoundaries: retrieve_interior_active_cells_map, ActiveCellsIBG, 
                                        active_linear_index_to_tuple
 
@@ -88,7 +90,8 @@ function compute_hydrostatic_free_surface_tendency_contributions!(model, kernel_
                      model.diffusivity_fields,
                      model.auxiliary_fields,
                      c_forcing,
-                     model.clock)
+                     model.clock,
+                     Val(compute_bgc_with_physics(model.timestepper)))
 
         for parameters in kernel_parameters
             launch!(arch, grid, parameters,
diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
index 6de795e2b6..264dea83f2 100644
--- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
+++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_tendency_kernel_functions.jl
@@ -117,7 +117,8 @@ where `c = C[tracer_index]`.
                                                           diffusivities,
                                                           auxiliary_fields,
                                                           forcing,
-                                                          clock) where tracer_index
+                                                          clock,
+                                                          val_compute_bgc) where tracer_index
 
     @inbounds c = tracers[tracer_index]
     model_fields = merge(hydrostatic_fields(velocities, free_surface, tracers), auxiliary_fields)
@@ -133,7 +134,7 @@ where `c = C[tracer_index]`.
     return ( - div_Uc(i, j, k, grid, advection, total_velocities, c)
              - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivities, val_tracer_index, c, clock, model_fields, buoyancy)
              - immersed_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, diffusivities, val_tracer_index, clock, model_fields)
-             + biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields)
+             + biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields, val_compute_bgc)
              + forcing(i, j, k, grid, clock, model_fields))
 end
 
diff --git a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl
index fb7e89ee64..1b6d32df26 100644
--- a/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl
+++ b/src/Models/NonhydrostaticModels/compute_nonhydrostatic_tendencies.jl
@@ -6,6 +6,8 @@ using Oceananigans.Models: complete_communication_and_compute_boundary!, interio
 using Oceananigans.ImmersedBoundaries: retrieve_interior_active_cells_map, ActiveCellsIBG, 
                                        active_linear_index_to_tuple
 
+using Oceananigans.TimeSteppers: compute_bgc_with_physics, timestepper_tendencies
+
 import Oceananigans.TimeSteppers: compute_tendencies!
 
 """
@@ -34,7 +36,7 @@ function compute_tendencies!(model::NonhydrostaticModel, callbacks)
 
     # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
     # boundaries of the domain
-    compute_boundary_tendency_contributions!(model.timestepper.Gⁿ,
+    compute_boundary_tendency_contributions!(timestepper_tendencies(model.timestepper),
                                              model.architecture,
                                              model.velocities,
                                              model.tracers,
@@ -53,7 +55,7 @@ end
 """ Store previous value of the source term and compute current source term. """
 function compute_interior_tendency_contributions!(model, kernel_parameters; active_cells_map = nothing)
 
-    tendencies           = model.timestepper.Gⁿ
+    tendencies           = timestepper_tendencies(model.timestepper)
     arch                 = model.architecture
     grid                 = model.grid
     advection            = model.advection
@@ -126,7 +128,8 @@ function compute_interior_tendency_contributions!(model, kernel_parameters; acti
                      start_tracer_kernel_args..., 
                      c_immersed_bc,
                      end_tracer_kernel_args...,
-                     forcing, clock)
+                     forcing, clock,
+                     Val(compute_bgc_with_physics(model.timestepper)))
 
         for parameters in kernel_parameters
             launch!(arch, grid, parameters, compute_Gc!, 
diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl
index cefb5eadc3..81744904c3 100644
--- a/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl
+++ b/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl
@@ -248,7 +248,8 @@ velocity components, tracer fields, and precalculated diffusivities where applic
                                  auxiliary_fields,
                                  diffusivities,
                                  forcing,
-                                 clock) where tracer_index
+                                 clock,
+                                 val_compute_bgc) where tracer_index
 
     @inbounds c = tracers[tracer_index]
     @inbounds background_fields_c = background_fields.tracers[tracer_index]
@@ -266,6 +267,6 @@ velocity components, tracer fields, and precalculated diffusivities where applic
              - div_Uc(i, j, k, grid, advection, velocities, background_fields_c)
              - ∇_dot_qᶜ(i, j, k, grid, closure, diffusivities, val_tracer_index, c, clock, model_fields, buoyancy)
              - immersed_∇_dot_qᶜ(i, j, k, grid, c, c_immersed_bc, closure, diffusivities, val_tracer_index, clock, model_fields)
-             + biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields)
+             + biogeochemical_transition(i, j, k, grid, biogeochemistry, val_tracer_name, clock, model_fields, val_compute_bgc)
              + forcing(i, j, k, grid, clock, model_fields))
 end
diff --git a/src/Models/ShallowWaterModels/rk3_substep_shallow_water_model.jl b/src/Models/ShallowWaterModels/rk3_substep_shallow_water_model.jl
index 880ae39163..2997cdb99c 100644
--- a/src/Models/ShallowWaterModels/rk3_substep_shallow_water_model.jl
+++ b/src/Models/ShallowWaterModels/rk3_substep_shallow_water_model.jl
@@ -14,14 +14,14 @@ function rk3_substep!(model::ShallowWaterModel, Δt, γⁿ, ζⁿ)
 
     substep_solution_kernel!(model.solution,
                              Δt, γⁿ, ζⁿ,
-                             model.timestepper.Gⁿ,
-                             model.timestepper.G⁻)
+                             timestepper_tendencies(timestepper),
+                             timestepper_previous_tendencies(timestepper))
 
 
     for i in 1:length(model.tracers)
         @inbounds c = model.tracers[i]
-        @inbounds Gcⁿ = model.timestepper.Gⁿ[i+3]
-        @inbounds Gc⁻ = model.timestepper.G⁻[i+3]
+        @inbounds Gcⁿ = timestepper_tendencies(timestepper)[i+3]
+        @inbounds Gc⁻ = timestepper_previous_tendencies(timestepper)[i+3]
 
         substep_tracer_kernel!(c, Δt, γⁿ, ζⁿ, Gcⁿ, Gc⁻)
     end
diff --git a/src/Models/interleave_communication_and_computation.jl b/src/Models/interleave_communication_and_computation.jl
index 286097e826..5d1e46fa21 100644
--- a/src/Models/interleave_communication_and_computation.jl
+++ b/src/Models/interleave_communication_and_computation.jl
@@ -6,6 +6,8 @@ using Oceananigans.DistributedComputations
 using Oceananigans.DistributedComputations: DistributedGrid
 using Oceananigans.DistributedComputations: synchronize_communication!, SynchronizedDistributed
 
+import Oceananigans: interior_tendency_kernel_parameters
+
 function complete_communication_and_compute_boundary!(model, ::DistributedGrid, arch)
 
     # Iterate over the fields to clear _ALL_ possible architectures
diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl
index 2abd97e7a0..036b8e204d 100644
--- a/src/Oceananigans.jl
+++ b/src/Oceananigans.jl
@@ -190,6 +190,8 @@ function fields end
 function prognostic_fields end
 function tracer_tendency_kernel_function end
 function boundary_conditions end
+function compute_biogeochemical_tendencies! end
+function interior_tendency_kernel_parameters end
 
 #####
 ##### Include all the submodules
diff --git a/src/TimeSteppers/TimeSteppers.jl b/src/TimeSteppers/TimeSteppers.jl
index 2955ca6999..ce791de1b9 100644
--- a/src/TimeSteppers/TimeSteppers.jl
+++ b/src/TimeSteppers/TimeSteppers.jl
@@ -55,9 +55,13 @@ step_lagrangian_particles!(model, Δt) = nothing
 reset!(timestepper) = nothing
 implicit_step!(field, ::Nothing, args...; kwargs...) = nothing
 
+compute_bgc_with_physics(timestepper) = true
+timestepper_tendencies(timestepper) = timestepper.Gⁿ
+timestepper_previous_tendencies(timestepper) = timestepper.G⁻
+
 include("clock.jl")
-include("store_tendencies.jl")
 include("quasi_adams_bashforth_2.jl")
 include("runge_kutta_3.jl")
-
+include("strange_splitting.jl")
+include("store_tendencies.jl")
 end # module
diff --git a/src/TimeSteppers/quasi_adams_bashforth_2.jl b/src/TimeSteppers/quasi_adams_bashforth_2.jl
index ec047f4a61..48a88ec3ab 100644
--- a/src/TimeSteppers/quasi_adams_bashforth_2.jl
+++ b/src/TimeSteppers/quasi_adams_bashforth_2.jl
@@ -73,16 +73,15 @@ The steps of the Quasi-Adams-Bashforth second-order (AB2) algorithm are:
 6. Update the model state.
 7. Compute tendencies for the next time step
 """
-function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt;
-                    callbacks=[], euler=false)
+time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; callbacks=[], euler=false) = 
+    time_step!(model.timestepper, model, Δt; callbacks, euler)
 
+function time_step!(ab2_timestepper::QuasiAdamsBashforth2TimeStepper, model, Δt; callbacks=[], euler=false)
     Δt == 0 && @warn "Δt == 0 may cause model blowup!"
 
     # Be paranoid and update state at iteration 0
     model.clock.iteration == 0 && update_state!(model, callbacks)
 
-    ab2_timestepper = model.timestepper
-
     # Change the default χ if necessary, which occurs if:
     #   * We detect that the time-step size has changed.
     #   * We detect that this is the "first" time-step, which means we
@@ -151,8 +150,8 @@ function ab2_step!(model, Δt)
     for (i, field) in enumerate(model_fields)
 
         step_field_kernel!(field, Δt, χ,
-                           model.timestepper.Gⁿ[i],
-                           model.timestepper.G⁻[i])
+                           timestepper_tendencies(timestepper)[i],
+                           timestepper_previous_tendencies(timestepper)[i])
 
         # TODO: function tracer_index(model, field_index) = field_index - 3, etc...
         tracer_index = Val(i - 3) # assumption
diff --git a/src/TimeSteppers/runge_kutta_3.jl b/src/TimeSteppers/runge_kutta_3.jl
index 0bc5c34608..7bc6830d69 100644
--- a/src/TimeSteppers/runge_kutta_3.jl
+++ b/src/TimeSteppers/runge_kutta_3.jl
@@ -1,4 +1,6 @@
 using Oceananigans.Architectures: architecture
+using Oceananigans: compute_biogeochemical_tendencies!
+
 using Oceananigans: fields
 
 """
@@ -78,23 +80,29 @@ The 3rd-order Runge-Kutta method takes three intermediate substep stages to
 achieve a single timestep. A pressure correction step is applied at each intermediate
 stage.
 """
-function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbacks=[])
+time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbacks=[]) =
+    time_step!(model.timestepper, model, Δt; callbacks)
+
+function time_step!(timestepper, model, Δt; callbacks=[])
     Δt == 0 && @warn "Δt == 0 may cause model blowup!"
 
     # Be paranoid and update state at iteration 0, in case run! is not used:
     model.clock.iteration == 0 && update_state!(model, callbacks; compute_tendencies = true)
 
-    γ¹ = model.timestepper.γ¹
-    γ² = model.timestepper.γ²
-    γ³ = model.timestepper.γ³
+    γ¹ = timestepper.γ¹
+    γ² = timestepper.γ²
+    γ³ = timestepper.γ³
 
-    ζ² = model.timestepper.ζ²
-    ζ³ = model.timestepper.ζ³
+    ζ² = timestepper.ζ²
+    ζ³ = timestepper.ζ³
 
     first_stage_Δt  = γ¹ * Δt
     second_stage_Δt = (γ² + ζ²) * Δt
     third_stage_Δt  = (γ³ + ζ³) * Δt
 
+    Gⁿ = timestepper_tendencies(timestepper)
+    G⁻ = timestepper_previous_tendencies(timestepper)
+    
     # Compute the next time step a priori to reduce floating point error accumulation
     tⁿ⁺¹ = next_time(model.clock, Δt)
 
@@ -102,7 +110,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac
     # First stage
     #
 
-    rk3_substep!(model, Δt, γ¹, nothing)
+    rk3_substep!(model, Δt, γ¹, nothing, Gⁿ, G⁻, timestepper.implicit_solver)
 
     tick!(model.clock, first_stage_Δt; stage=true)
     model.clock.last_stage_Δt = first_stage_Δt
@@ -118,7 +126,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac
     # Second stage
     #
 
-    rk3_substep!(model, Δt, γ², ζ²)
+    rk3_substep!(model, Δt, γ², ζ², Gⁿ, G⁻, timestepper.implicit_solver)
 
     tick!(model.clock, second_stage_Δt; stage=true)
     model.clock.last_stage_Δt = second_stage_Δt
@@ -134,7 +142,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac
     # Third stage
     #
     
-    rk3_substep!(model, Δt, γ³, ζ³)
+    rk3_substep!(model, Δt, γ³, ζ³, Gⁿ, G⁻, timestepper.implicit_solver)
 
     # This adjustment of the final time-step reduces the accumulation of
     # round-off error when Δt is added to model.clock.time. Note that we still use 
@@ -154,6 +162,74 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac
     return nothing
 end
 
+
+function time_step_biogeochemistry!(timestepper, model, Δt)
+    Δt == 0 && @warn "Δt == 0 may cause model blowup!"
+
+    # Be paranoid and update state at iteration 0, in case run! is not used:
+    model.clock.iteration == 0 && compute_biogeochemical_tendencies!(model, timestepper_tendencies(timestepper))
+
+    γ¹ = timestepper.γ¹
+    γ² = timestepper.γ²
+    γ³ = timestepper.γ³
+
+    ζ² = timestepper.ζ²
+    ζ³ = timestepper.ζ³
+
+    Gⁿ = timestepper_tendencies(timestepper)
+    G⁻ = timestepper_previous_tendencies(timestepper)
+
+    first_stage_Δt  = γ¹ * Δt
+    second_stage_Δt = (γ² + ζ²) * Δt
+    third_stage_Δt  = (γ³ + ζ³) * Δt
+
+    # Compute the next time step a priori to reduce floating point error accumulation
+    tⁿ⁺¹ = next_time(model.clock, Δt)
+
+    #
+    # First stage
+    #
+
+    rk3_substep!(model, Δt, γ¹, nothing, Gⁿ, G⁻, nothing)
+
+    tick!(model.clock, first_stage_Δt; stage=true)
+    model.clock.last_stage_Δt = first_stage_Δt
+
+    store_biogeochemical_tendencies!(model)
+    compute_biogeochemical_tendencies!(model, timestepper_tendencies(timestepper))
+
+    #
+    # Second stage
+    #
+
+    rk3_substep!(model, Δt, γ², ζ², Gⁿ, G⁻, nothing)
+
+    tick!(model.clock, second_stage_Δt; stage=true)
+    model.clock.last_stage_Δt = second_stage_Δt
+
+    store_biogeochemical_tendencies!(model)
+    compute_biogeochemical_tendencies!(model, timestepper_tendencies(timestepper))
+
+    #
+    # Third stage
+    #
+    
+    rk3_substep!(model, Δt, γ³, ζ³, Gⁿ, G⁻, nothing)
+
+    # This adjustment of the final time-step reduces the accumulation of
+    # round-off error when Δt is added to model.clock.time. Note that we still use 
+    # third_stage_Δt for the substep, pressure correction, and Lagrangian particles step.
+    corrected_third_stage_Δt = tⁿ⁺¹ - model.clock.time
+
+    tick!(model.clock, third_stage_Δt)
+    model.clock.last_stage_Δt = corrected_third_stage_Δt
+    model.clock.last_Δt = Δt
+
+    compute_biogeochemical_tendencies!(model, timestepper_tendencies(timestepper))
+
+    return nothing
+end
+
 #####
 ##### Time stepping in each substep
 #####
@@ -161,22 +237,20 @@ end
 stage_Δt(Δt, γⁿ, ζⁿ) = Δt * (γⁿ + ζⁿ)
 stage_Δt(Δt, γⁿ, ::Nothing) = Δt * γⁿ
 
-function rk3_substep!(model, Δt, γⁿ, ζⁿ)
+function rk3_substep!(model, Δt, γⁿ, ζⁿ, Gⁿ, G⁻, implicit_solver)
 
     workgroup, worksize = work_layout(model.grid, :xyz)
     substep_field_kernel! = rk3_substep_field!(device(architecture(model)), workgroup, worksize)
     model_fields = prognostic_fields(model)
 
     for (i, field) in enumerate(model_fields)
-        substep_field_kernel!(field, Δt, γⁿ, ζⁿ,
-                              model.timestepper.Gⁿ[i],
-                              model.timestepper.G⁻[i])
+        substep_field_kernel!(field, Δt, γⁿ, ζⁿ, Gⁿ[i], G⁻[i])
 
         # TODO: function tracer_index(model, field_index) = field_index - 3, etc...
         tracer_index = Val(i - 3) # assumption
 
         implicit_step!(field,
-                       model.timestepper.implicit_solver,
+                       implicit_solver,
                        model.closure,
                        model.diffusivity_fields,
                        tracer_index,
diff --git a/src/TimeSteppers/store_tendencies.jl b/src/TimeSteppers/store_tendencies.jl
index 2eae58dc2c..7e0d0a451e 100644
--- a/src/TimeSteppers/store_tendencies.jl
+++ b/src/TimeSteppers/store_tendencies.jl
@@ -14,8 +14,20 @@ function store_tendencies!(model)
 
     for field_name in keys(model_fields)
         launch!(model.architecture, model.grid, :xyz, store_field_tendencies!,
-                model.timestepper.G⁻[field_name],
-                model.timestepper.Gⁿ[field_name])
+                timestepper_previous_tendencies(model.timestepper)[field_name],
+                timestepper_tendencies(model.timestepper)[field_name])
+    end
+
+    return nothing
+end
+
+function store_biogeochemical_tendencies!(model::AbstractModel{<:StrangeSplittingTimeStepper})
+    model_fields = prognostic_fields(model)
+
+    for field_name in keys(model_fields)
+        launch!(model.architecture, model.grid, :xyz, store_field_tendencies!,
+                model.timestepper.biogeochemistry.G⁻[field_name],
+                model.timestepper.biogeochemistry.Gⁿ[field_name])
     end
 
     return nothing
diff --git a/src/TimeSteppers/strange_splitting.jl b/src/TimeSteppers/strange_splitting.jl
new file mode 100644
index 0000000000..7f9343021a
--- /dev/null
+++ b/src/TimeSteppers/strange_splitting.jl
@@ -0,0 +1,46 @@
+using Oceananigans.Fields: FunctionField, location
+using Oceananigans.Utils: @apply_regionally, apply_regionally!
+
+"""
+    StrangeSplittingTimeStepper
+
+Strange splitting 
+"""
+struct StrangeSplittingTimeStepper{PT, BT} <: AbstractTimeStepper
+            physics :: PT
+    biogeochemistry :: BT
+
+    function StrangeSplittingTimeStepper(grid, tracers;
+                                         implicit_solver = nothing,
+                                         physics::PT = RungeKutta3TimeStepper(grid, tracers; implicit_solver),
+                                         biogeochemistry::BT = RungeKutta3TimeStepper(grid, tracers)) where {PT, BT}
+        isa(biogeochemistry, RungeKutta3TimeStepper) || @warn "RK3 is the only biogeochemistry timestepper currently supported"
+
+        return new{PT, BT}(physics, biogeochemistry)
+    end
+end
+
+reset!(ts::StrangeSplittingTimeStepper) = (reset!(ts.physics); reset!(ts.biogeochemistry))
+
+"""
+    time_step!(model::AbstractModel{<:StrangeSplittingTimeStepper}, Δt;)
+
+"""
+function time_step!(model::AbstractModel{<:StrangeSplittingTimeStepper}, Δt; callbacks=[])
+    timesteppers = model.timestepper
+
+    time_step_biogeochemistry!(timesteppers.biogeochemistry, model, Δt/2)
+
+    time_step!(timesteppers.physics, model, Δt; callbacks)
+
+    time_step_biogeochemistry!(timesteppers.biogeochemistry, model, Δt/2)
+    
+    return nothing
+end
+
+compute_bgc_with_physics(::StrangeSplittingTimeStepper) = false
+
+timestepper_tendencies(timestepper::StrangeSplittingTimeStepper) = timestepper.physics.Gⁿ
+timestepper_previous_tendencies(timestepper::StrangeSplittingTimeStepper) = timestepper.physics.G⁻
+
+# TODO: implement step bgc tracers... - does this need to include surface fluxs? Probably, is that handelled by fill halo still?
\ No newline at end of file