diff --git a/Project.toml b/Project.toml
index df260ad..df8b4e0 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "Oceanostics"
 uuid = "d0ccf422-c8fb-49b5-a76d-74acdde946ac"
-version = "0.16.5"
+version = "0.16.6"
 authors = ["tomchor <tomaschor@gmail.com>"]
 
 [deps]
@@ -11,7 +11,7 @@ SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
 
 [compat]
 DocStringExtensions = "0.9"
-Oceananigans = "0.101"
+Oceananigans = "0.102"
 SeawaterPolynomials = "0.3"
 julia = "1.9"
 
diff --git a/src/FlowDiagnostics.jl b/src/FlowDiagnostics.jl
index 7272f7f..df5849b 100755
--- a/src/FlowDiagnostics.jl
+++ b/src/FlowDiagnostics.jl
@@ -18,7 +18,7 @@ using Oceananigans.BuoyancyFormulations: get_temperature_and_salinity, SeawaterB
 using Oceananigans.Operators
 using Oceananigans.AbstractOperations
 using Oceananigans.AbstractOperations: KernelFunctionOperation
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using Oceananigans.Grids: AbstractGrid, Center, Face, NegativeZDirection, ZDirection, znode, bottommost_active_node
 
 using SeawaterPolynomials: ρ′, BoussinesqEquationOfState
@@ -79,7 +79,7 @@ where `z` is the true vertical direction (ie anti-parallel to gravity).
 """
 function RichardsonNumber(model; loc = (Center, Center, Face))
     validate_location(loc, "RichardsonNumber", (Center, Center, Face))
-    return RichardsonNumber(model, model.velocities..., buoyancy(model); loc)
+    return RichardsonNumber(model, model.velocities..., buoyancy_operation(model); loc)
 end
 
 function RichardsonNumber(model, u, v, w, b; loc = (Center, Center, Face))
diff --git a/src/TracerEquation.jl b/src/TracerEquation.jl
index 980a0e2..84b5d04 100755
--- a/src/TracerEquation.jl
+++ b/src/TracerEquation.jl
@@ -194,7 +194,7 @@ julia> model = NonhydrostaticModel(; grid, tracers=:a);
 julia> FORC = TracerEquation.Forcing(model, :a)
 KernelFunctionOperation at (Center, Center, Center)
 ├── grid: 4×4×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
-├── kernel_function: zeroforcing (generic function with 1 method)
+├── kernel_function: Returns (generic function with 1 method)
 └── arguments: ("Clock", "NamedTuple")
 ```
 """
diff --git a/src/TracerVarianceEquation.jl b/src/TracerVarianceEquation.jl
index e7bb6ba..1f3357d 100644
--- a/src/TracerVarianceEquation.jl
+++ b/src/TracerVarianceEquation.jl
@@ -62,7 +62,7 @@ julia> χ = TracerVarianceEquation.TracerVarianceTendency(model, :b)
 KernelFunctionOperation at (Center, Center, Center)
 ├── grid: 1×1×4 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo
 ├── kernel_function: c∂ₜcᶜᶜᶜ (generic function with 1 method)
-└── arguments: ("Val", "Val", "Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Clock", "typeof(Oceananigans.Forcings.zeroforcing)")
+└── arguments: ("Val", "Val", "Centered", "Nothing", "Nothing", "Nothing", "Nothing", "Oceananigans.Models.NonhydrostaticModels.BackgroundFields", "NamedTuple", "NamedTuple", "NamedTuple", "Nothing", "Clock", "Returns")
 ```
 """
 function TracerVarianceTendency(model::NonhydrostaticModel, tracer_name; location = (Center, Center, Center))
diff --git a/test/test_active_tracer_diagnostics.jl b/test/test_active_tracer_diagnostics.jl
index 3f7f215..de5542f 100644
--- a/test/test_active_tracer_diagnostics.jl
+++ b/test/test_active_tracer_diagnostics.jl
@@ -3,7 +3,7 @@ using CUDA: @allowscalar
 using Oceananigans: fill_halo_regions!
 using Oceananigans.AbstractOperations: AbstractOperation
 using Oceananigans.Grids: znode
-using Oceananigans.BuoyancyFormulations: buoyancy
+using Oceananigans.Models: buoyancy_operation
 using SeawaterPolynomials: BoussinesqEquationOfState
 using SeawaterPolynomials.SecondOrderSeawaterPolynomials: LinearRoquetSeawaterPolynomial
 
@@ -47,7 +47,7 @@ function test_buoyancy_diagnostics(model)
         @test Ri_field isa Field
         @test interior(Ri_field, 3, 3, 3)[1] ≈ N² / S^2
 
-        b = buoyancy(model)
+        b = buoyancy_operation(model)
         Ri = RichardsonNumber(model, u, v, w, b)
         @test Ri isa RichardsonNumber
         Ri_field = Field(Ri)