diff --git a/Project.toml b/Project.toml
index 441db7d7f..23b0f637a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
-version = "0.8.7"
+version = "0.8.8"
 authors = ["Climate Modeling Alliance and contributors"]
 
 [deps]
@@ -60,7 +60,7 @@ KernelAbstractions = "0.9"
 MPI = "0.20"
 MeshArrays = "0.3, 0.4"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.101"
+Oceananigans = "0.102"
 OffsetArrays = "1.14"
 PrecompileTools = "1"
 Reactant = "0.2.45"
diff --git a/examples/ecco_inspect_temperature_salinity.jl b/examples/ecco_inspect_temperature_salinity.jl
index df1e2210f..026cbbdd2 100644
--- a/examples/ecco_inspect_temperature_salinity.jl
+++ b/examples/ecco_inspect_temperature_salinity.jl
@@ -31,12 +31,12 @@ T = CenterField(grid)
 S = CenterField(grid)
 
 using SeawaterPolynomials: TEOS10EquationOfState
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 
 equation_of_state = TEOS10EquationOfState()
 sb = SeawaterBuoyancy(; equation_of_state)
 tracers = (T=T, S=S)
-b = Field(buoyancy(sb, grid, tracers))
+b = Field(buoyancy_operation(sb, grid, tracers))
 
 start_date = DateTime(1993, 1, 1)
 end_date   = DateTime(1999, 1, 1)
diff --git a/examples/ecco_mixed_layer_depth.jl b/examples/ecco_mixed_layer_depth.jl
index 78217ba37..cead5fc61 100644
--- a/examples/ecco_mixed_layer_depth.jl
+++ b/examples/ecco_mixed_layer_depth.jl
@@ -6,7 +6,7 @@ using Printf
 using Dates
 
 using SeawaterPolynomials: TEOS10EquationOfState
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 
 arch = CPU()
 Nx = 360
@@ -44,7 +44,7 @@ Nt = length(ht)
 for n = 1:Nt-1
     local tracers
     tracers = (T=Tt[n], S=St[n])
-    h.operand.buoyancy_perturbation = buoyancy(sb, grid, tracers)
+    h.operand.buoyancy_perturbation = buoyancy_operation(sb, grid, tracers)
     @show n
     @time compute!(h)
     parent(ht[n]) .= parent(h)
diff --git a/examples/single_column_os_papa_simulation.jl b/examples/single_column_os_papa_simulation.jl
index 5c59b4a82..6b92168a3 100644
--- a/examples/single_column_os_papa_simulation.jl
+++ b/examples/single_column_os_papa_simulation.jl
@@ -18,7 +18,6 @@ using ClimaOcean
 using Oceananigans
 using Oceananigans.Units
 using Oceananigans.BuoyancyFormulations: buoyancy_frequency
-using Oceananigans.Units: Time
 using Dates
 using Printf
 
diff --git a/experiments/one_degree_calibration/distributed_perfect_one_degree_model_calibration.jl b/experiments/one_degree_calibration/distributed_perfect_one_degree_model_calibration.jl
index e23b6b532..18dc5f08e 100644
--- a/experiments/one_degree_calibration/distributed_perfect_one_degree_model_calibration.jl
+++ b/experiments/one_degree_calibration/distributed_perfect_one_degree_model_calibration.jl
@@ -1,9 +1,10 @@
 using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Utils: WallTimeInterval
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
 using ClimaOcean.NearGlobalSimulations: one_degree_near_global_simulation
+
 using ParameterEstimocean
 using ParameterEstimocean.Utils: map_gpus_to_ranks!
 using ParameterEstimocean.Observations: FieldTimeSeriesCollector
diff --git a/experiments/one_degree_calibration/gm_one_degree_model_calibration.jl b/experiments/one_degree_calibration/gm_one_degree_model_calibration.jl
index 1b5c5bb83..674fe9def 100644
--- a/experiments/one_degree_calibration/gm_one_degree_model_calibration.jl
+++ b/experiments/one_degree_calibration/gm_one_degree_model_calibration.jl
@@ -2,7 +2,7 @@ using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Grids: on_architecture
 using Oceananigans.Utils: WallTimeInterval
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
 using ClimaOcean.NearGlobalSimulations: one_degree_near_global_simulation
 using ParameterEstimocean
diff --git a/experiments/one_degree_calibration/one_degree_near_global_simulation.jl b/experiments/one_degree_calibration/one_degree_near_global_simulation.jl
index d37d50f05..86c09e299 100644
--- a/experiments/one_degree_calibration/one_degree_near_global_simulation.jl
+++ b/experiments/one_degree_calibration/one_degree_near_global_simulation.jl
@@ -3,7 +3,7 @@ using ClimaOcean.NearGlobalSimulations: one_degree_near_global_simulation
 using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Utils: WallTimeInterval
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
 using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
     MixingLength, TurbulentKineticEnergyEquation, CATKEVerticalDiffusivity
diff --git a/src/DataWrangling/metadata.jl b/src/DataWrangling/metadata.jl
index cca37a744..3b45eed05 100644
--- a/src/DataWrangling/metadata.jl
+++ b/src/DataWrangling/metadata.jl
@@ -1,7 +1,7 @@
 using CFTime
 using Dates
 using Base: @propagate_inbounds
-import Oceananigans.Utils: prettysummary
+import Oceananigans.Grids: prettysummary
 
 struct BoundingBox{X, Y, Z}
     longitude :: X
diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl
index 489efb1d0..ad5db51a6 100644
--- a/src/DataWrangling/restoring.jl
+++ b/src/DataWrangling/restoring.jl
@@ -2,7 +2,7 @@ using Oceananigans: location
 using Oceananigans.Grids: node
 using Oceananigans.Fields: interpolate, instantiated_location
 using Oceananigans.OutputReaders: Cyclical
-using Oceananigans.Utils: Time
+using Oceananigans.Units: Time
 using Oceananigans.Architectures: AbstractArchitecture
 
 using JLD2
diff --git a/src/Diagnostics/Diagnostics.jl b/src/Diagnostics/Diagnostics.jl
index 7f3931f12..76e21b8c0 100644
--- a/src/Diagnostics/Diagnostics.jl
+++ b/src/Diagnostics/Diagnostics.jl
@@ -4,16 +4,15 @@ export MixedLayerDepthField, MixedLayerDepthOperand
 
 using Oceananigans
 using Oceananigans.Architectures: architecture
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Grids: new_data, inactive_cell, znode
 using Oceananigans.BoundaryConditions: FieldBoundaryConditions, fill_halo_regions!
 using Oceananigans.Fields: FieldStatus
 using Oceananigans.Utils: launch!
+using KernelAbstractions: @index, @kernel
 
 import Oceananigans.Fields: compute!
 
-using KernelAbstractions: @index, @kernel
-
 include("mixed_layer_depth.jl")
 
 end # module
diff --git a/src/Diagnostics/mixed_layer_depth.jl b/src/Diagnostics/mixed_layer_depth.jl
index 3b8b81a65..0b986312b 100644
--- a/src/Diagnostics/mixed_layer_depth.jl
+++ b/src/Diagnostics/mixed_layer_depth.jl
@@ -8,7 +8,7 @@ end
 Base.summary(mldo::MixedLayerDepthOperand) = "MixedLayerDepthOperand"
 
 function MixedLayerDepthOperand(bm, grid, tracers; difference_criterion=1e-4)
-    buoyancy_perturbation = buoyancy(bm, grid, tracers)
+    buoyancy_perturbation = buoyancy_operation(bm, grid, tracers)
     difference_criterion = convert(eltype(grid), difference_criterion)
     return MixedLayerDepthOperand(buoyancy_perturbation, difference_criterion)
 end
@@ -16,7 +16,7 @@ end
 const MixedLayerDepthField = Field{<:Any, <:Any, <:Any, <:MixedLayerDepthOperand}
 
 """
-    MixedLayerDepthField(bm, grid, tracers; difference_criterion=1e-4)
+    MixedLayerDepthField(bm, grid, tracers; difference_criterion=3e-5)
 
 """
 function MixedLayerDepthField(bm, grid, tracers; difference_criterion=3e-5)
@@ -92,4 +92,3 @@ const f = Face()
 
     @inbounds mld[i, j, 1] = h
 end
-
diff --git a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
index 215a9bcaa..32da8cd9d 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl
@@ -17,10 +17,11 @@ using ..OceanSeaIceModels.PrescribedAtmospheres:
 using ClimaSeaIce: SeaIceModel
 
 using Oceananigans: HydrostaticFreeSurfaceModel, architecture
+using Oceananigans.Units: Time
 using Oceananigans.Grids: inactive_node, node, topology
 using Oceananigans.BoundaryConditions: fill_halo_regions!
 using Oceananigans.Fields: ConstantField, interpolate, FractionalIndices
-using Oceananigans.Utils: launch!, Time, KernelParameters
+using Oceananigans.Utils: launch!, KernelParameters
 using Oceananigans.Operators: ℑxᶜᵃᵃ, ℑyᵃᶜᵃ, ℑxᶠᵃᵃ, ℑyᵃᶠᵃ
 
 using KernelAbstractions: @kernel, @index
diff --git a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl
index 469e5c827..fc657b7f4 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl
@@ -1,5 +1,4 @@
-using Oceananigans.Utils: prettysummary
-using Oceananigans.Grids: AbstractGrid
+using Oceananigans.Grids: AbstractGrid, prettysummary
 
 using Adapt
 using Printf
diff --git a/src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl b/src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl
index 9607ded5a..2007f6a2f 100644
--- a/src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl
+++ b/src/OceanSeaIceModels/InterfaceComputations/tabulated_albedo.jl
@@ -1,6 +1,6 @@
 using Oceananigans.Fields: interpolator
 using Oceananigans.Grids: on_architecture
-using Oceananigans.Utils: Time
+using Oceananigans.Units: Time
 using Base
 
 using ClimaOcean.OceanSeaIceModels:
diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl
index a5226076b..d024dd3db 100644
--- a/src/OceanSeaIceModels/OceanSeaIceModels.jl
+++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl
@@ -14,7 +14,8 @@ using SeawaterPolynomials
 
 using Oceananigans
 using Oceananigans.Operators
-using Oceananigans.Utils: launch!, Time, KernelParameters
+using Oceananigans.Utils: launch!, KernelParameters
+using Oceananigans.Units: Time
 using Oceananigans.Architectures: architecture
 using Oceananigans.BoundaryConditions: fill_halo_regions!, BoundaryCondition
 using Oceananigans.Grids: architecture
diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl
index 49d72d32d..46a9c5de0 100644
--- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl
+++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl
@@ -5,7 +5,8 @@ using Oceananigans.Fields: Center
 using Oceananigans.Grids: grid_name
 using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series
 using Oceananigans.TimeSteppers: Clock, tick!
-using Oceananigans.Utils: prettysummary, Time
+using Oceananigans.Grids: prettysummary
+using Oceananigans.Units: Time
 
 using Adapt
 using Thermodynamics.Parameters: AbstractThermodynamicsParameters
diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl
index ed95c4f14..b532fa00a 100644
--- a/test/runtests_setup.jl
+++ b/test/runtests_setup.jl
@@ -198,16 +198,11 @@ function test_timestepping_with_dataset_restoring(arch, dataset, dates, inpainti
                                  z = (-200, 0),
                                  halo = (6, 6, 6))
 
-    # Dynamically create name of forcing based on dataset field name
-    forcing = NamedTuple{
-                (fldnames)
-            }(
-                ntuple(i->DatasetRestoring(
-                                           Metadata(varnames[i]; dates, dataset), 
-                                           arch; inpainting, rate=1/1000
-                                          ), length(varnames))
-             )
-
+    # Force only the last tracer. 
+    # Forcing more than one variable leads to parameter space errors
+    metadata = Metadata(varnames[end]; dates, dataset)
+    restoring = DatasetRestoring(metadata, arch; inpainting, rate=1/1000)
+    forcing = NamedTuple{tuple(fldnames[end])}(tuple(restoring))
     ocean = ocean_simulation(grid; tracers=fldnames, forcing, verbose=false)
     
     @test begin
@@ -234,50 +229,47 @@ function test_cycling_dataset_restoring(arch, dataset, dates, inpainting;
     start_date = dates[1]
     end_date = dates[end]
 
-    metadata1 = Metadata(varnames[end]; dates, dataset)
-    metadata2 = Metadata(varnames[end]; start_date, end_date, dataset)
+    metadata = Metadata(varnames[end]; dates, dataset)
 
-    for metadata in (metadata1, metadata2)
-        # Dynamically create name of forcing based on dataset field name
-        # Dynamically create name of forcing based on dataset field name
-        forcing = NamedTuple{
-                (fldnames[end],)
-            }(
-                (DatasetRestoring(metadata, arch;  time_indices_in_memory, inpainting, rate=1/1000),)
-             )
+    # Dynamically create name of forcing based on dataset field name
+    # Dynamically create name of forcing based on dataset field name
+    forcing = NamedTuple{
+            (fldnames[end],)
+        }(
+            (DatasetRestoring(metadata, arch;  time_indices_in_memory, inpainting, rate=1/1000),)
+            )
 
-        times = native_times(forcing[1].field_time_series.backend.metadata)
-        ocean = ocean_simulation(grid, tracers=fldnames, forcing=forcing)
+    times = native_times(forcing[1].field_time_series.backend.metadata)
+    ocean = ocean_simulation(grid, tracers=fldnames, forcing=forcing)
 
-        # start a bit after time_index
-        time_index = 3
-        time_interval = dataset isa ECCO2Daily ? Units.hour : Units.day
-        ocean.model.clock.time = times[time_index] + 2 * time_interval
-        update_state!(ocean.model)
+    # start a bit after time_index
+    time_index = 3
+    time_interval = dataset isa ECCO2Daily ? Units.hour : Units.day
+    ocean.model.clock.time = times[time_index] + 2 * time_interval
+    update_state!(ocean.model)
 
-        @test time_indices(forcing[1].field_time_series) ==
-            Tuple(range(time_index, length=time_indices_in_memory))
+    @test time_indices(forcing[1].field_time_series) ==
+        Tuple(range(time_index, length=time_indices_in_memory))
 
-        @test forcing[1].field_time_series.backend.start == time_index
+    @test forcing[1].field_time_series.backend.start == time_index
 
-        # Compile
-        time_step!(ocean)
+    # Compile
+    time_step!(ocean)
 
-        # Try stepping out of the dataset bounds
-        # start a bit after last time_index
-        ocean.model.clock.time = last(times) + 2 * time_interval
+    # Try stepping out of the dataset bounds
+    # start a bit after last time_index
+    ocean.model.clock.time = last(times) + 2 * time_interval
 
-        update_state!(ocean.model)
+    update_state!(ocean.model)
 
-        @test begin
-            time_step!(ocean)
-            true
-        end
-
-        # The backend has cycled to the end
-        @test time_indices(forcing[1].field_time_series) ==
-            mod1.(Tuple(range(length(times), length=time_indices_in_memory)), length(times))
+    @test begin
+        time_step!(ocean)
+        true
     end
+
+    # The backend has cycled to the end
+    @test time_indices(forcing[1].field_time_series) ==
+        mod1.(Tuple(range(length(times), length=time_indices_in_memory)), length(times))
 end
 
 function test_inpainting_algorithm(arch, dataset, start_date, inpainting; 
diff --git a/test/test_diagnostics.jl b/test/test_diagnostics.jl
index df788ec3e..62708b6a6 100644
--- a/test/test_diagnostics.jl
+++ b/test/test_diagnostics.jl
@@ -1,8 +1,8 @@
 include("runtests_setup.jl")
 
 using SeawaterPolynomials: TEOS10EquationOfState
-using Oceananigans.BuoyancyFormulations: buoyancy
 using Oceananigans: location
+using Oceananigans.Models: buoyancy_operation
 using ClimaOcean.Diagnostics: MixedLayerDepthField, MixedLayerDepthOperand
 
 for arch in test_architectures, dataset in (ECCO4Monthly(),)
@@ -49,7 +49,7 @@ for arch in test_architectures, dataset in (ECCO4Monthly(),)
         end
 
         tracers = (T=Tt[2], S=St[2])
-        h.operand.buoyancy_perturbation = buoyancy(sb, grid, tracers)
+        h.operand.buoyancy_perturbation = buoyancy_operation(sb, grid, tracers)
         compute!(h)
         if dataset isa ECCO4Monthly
             @test @allowscalar h[1, 1, 1] ≈ 9.2957298 # m
diff --git a/test/test_ecco2_monthly.jl b/test/test_ecco2_monthly.jl
index a74d9d743..ad4fbeb2f 100644
--- a/test/test_ecco2_monthly.jl
+++ b/test/test_ecco2_monthly.jl
@@ -59,7 +59,7 @@ for arch in test_architectures, dataset in test_ecco_datasets
 
             @testset "Field utilities" begin
                 test_ocean_metadata_utilities(arch, dataset, dates, inpainting,
-                                            varnames=test_names[dataset])
+                                              varnames=test_names[dataset])
             end
 
             @testset "DatasetRestoring with LinearlyTaperedPolarMask" begin