diff --git a/Project.toml b/Project.toml
index 90404515b1..72386a5104 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Oceananigans"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.100.7"
+version = "0.100.8"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl
index d24cdd6b7c..5c0d09db04 100644
--- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl
+++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl
@@ -1,3 +1,5 @@
+using Oceananigans.Utils: time_difference_seconds
+
 struct CATKEVerticalDiffusivity{TD, CL, FT, DT, TKE} <: AbstractScalarDiffusivity{TD, VerticalFormulation, 2}
     mixing_length :: CL
     turbulent_kinetic_energy_equation :: TKE
@@ -223,7 +225,7 @@ end
 @inline diffusivity_location(::FlavorOfCATKE) = (c, c, f)
 
 function update_previous_compute_time!(diffusivities, model)
-    Δt = model.clock.time - diffusivities.previous_compute_time[]
+    Δt = time_difference_seconds(model.clock.time, diffusivities.previous_compute_time[])
     diffusivities.previous_compute_time[] = model.clock.time
     return Δt
 end
diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl
index 684bbb87b6..df281373c3 100644
--- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl
+++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl
@@ -55,6 +55,8 @@ struct VariableStabilityFunctions{FT} <: AbstractConstantSchmidtStabilityFunctio
     𝕊u₀ :: FT
 end
 
+VariableStabilityFunctions{FT}(; kw...) where FT = VariableStabilityFunctions(FT; kw...)
+
 function VariableStabilityFunctions(FT=Oceananigans.defaults.FloatType;
                                     Cσe = 1.0,
                                     Cσϵ = 1.2,
@@ -75,7 +77,7 @@ function VariableStabilityFunctions(FT=Oceananigans.defaults.FloatType;
     if isnothing(𝕊u₀)
         # Compute 𝕊u₀ for the logarithmic boundary layer where production
         # balances dissipation. For more information see the discussion
-        # surrounding equation (13) in Umlauf and Burchard (2003).
+        # surrounding equation (13) by Umlauf and Burchard (2003).
         a = Cd₅ - Cu₂
         b = Cd₂ - Cu₀
         c = Cd₀
diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl
index e4bdb42e5c..e44be9ebc8 100644
--- a/test/test_time_stepping.jl
+++ b/test/test_time_stepping.jl
@@ -309,22 +309,33 @@ timesteppers = (:QuasiAdamsBashforth2, :RungeKutta3)
     @info "Testing time stepping..."
 
     for arch in archs, FT in float_types
-        @testset "Time stepping with DateTimes [$(typeof(arch)), $FT]" begin
-            @info "  Testing time stepping with datetime clocks [$(typeof(arch)), $FT]"
+        A = typeof(arch)
+        @testset "Time stepping with DateTimes [$A, $FT]" begin
+            @info "  Testing NonhydrostaticModel time stepping with datetime clocks [$A, $FT]"
 
             grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
             clock = Clock(time=DateTime(2020))
-            model = NonhydrostaticModel(; grid, clock, timestepper=:QuasiAdamsBashforth2)
+            model = NonhydrostaticModel(; grid, clock)
 
             time_step!(model, 7.883)
             @test model.clock.time == DateTime("2020-01-01T00:00:07.883")
 
-            model = NonhydrostaticModel(grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1)),
-                                        timestepper = :QuasiAdamsBashforth2,
-                                        clock = Clock(time=TimeDate(2020)))
-
+            clock = Clock(time=TimeDate(2020))
+            model = NonhydrostaticModel(; grid, clock)
             time_step!(model, 123e-9)  # 123 nanoseconds
             @test model.clock.time == TimeDate("2020-01-01T00:00:00.000000123")
+
+            # Test HydrostaticFreeSurfaceModel
+            for closure in (nothing, CATKEVerticalDiffusivity(FT), TKEDissipationVerticalDiffusivity(FT))
+                @info "  Testing HydrostaticFreeSurfaceModel time stepping with datetime clocks [$A, $FT, $(typeof(closure))]"
+                tracers = (:b, :c, :e, :ϵ)
+                clock = Clock(time=DateTime(2020, 1, 1))
+                grid = RectilinearGrid(arch, FT; size=(2, 2, 2), extent=(1, 1, 1))
+                model = HydrostaticFreeSurfaceModel(; grid, clock, closure, tracers,
+                                                    buoyancy = BuoyancyTracer())
+                time_step!(model, 1)
+                @test model.clock.time == DateTime("2020-01-01T00:00:01")
+            end
         end
     end