diff --git a/docs/src/model_setup/buoyancy_and_equation_of_state.md b/docs/src/model_setup/buoyancy_and_equation_of_state.md
index bc2e589ee9..8329573b8d 100644
--- a/docs/src/model_setup/buoyancy_and_equation_of_state.md
+++ b/docs/src/model_setup/buoyancy_and_equation_of_state.md
@@ -156,9 +156,9 @@ we might alternatively specify
 julia> buoyancy = SeawaterBuoyancy(gravitational_acceleration=1.3)
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 1.3
-└── equation of state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078)
+└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078)
 
-julia> model = NonhydrostaticModel(; grid, buoyancy=buoyancy, tracers=(:T, :S))
+julia> model = NonhydrostaticModel(; grid, buoyancy, tracers=(:T, :S))
 NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
 ├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
 ├── timestepper: QuasiAdamsBashforth2TimeStepper
@@ -179,7 +179,7 @@ To specify the thermal expansion and haline contraction coefficients
 julia> buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-3, haline_contraction=5e-4))
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
-└── equation of state: LinearEquationOfState(thermal_expansion=0.002, haline_contraction=0.0005)
+└── equation_of_state: LinearEquationOfState(thermal_expansion=0.002, haline_contraction=0.0005)
 ```
 
 ### Idealized nonlinear equations of state
@@ -202,7 +202,7 @@ julia> eos.seawater_polynomial # the density anomaly
 julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
-└── equation of state: BoussinesqEquationOfState{Float64}
+└── equation_of_state: BoussinesqEquationOfState{Float64}
 ```
 
 ### TEOS-10 equation of state
@@ -231,7 +231,7 @@ julia> θ = 45; # degrees
 
 julia> g̃ = (0, sind(θ), cosd(θ));
 
-julia> model = NonhydrostaticModel(; grid, 
+julia> model = NonhydrostaticModel(; grid,
                                    buoyancy=Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=g̃), 
                                    tracers=:b)
 NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
diff --git a/src/BuoyancyModels/seawater_buoyancy.jl b/src/BuoyancyModels/seawater_buoyancy.jl
index 88bd193a5f..ba7fa68390 100644
--- a/src/BuoyancyModels/seawater_buoyancy.jl
+++ b/src/BuoyancyModels/seawater_buoyancy.jl
@@ -36,7 +36,7 @@ function Base.show(io::IO, b::SeawaterBuoyancy{FT}) where FT
         print(io, "├── constant_salinity: ", b.constant_salinity, "\n")
     end
         
-    print(io, "└── equation of state: ", summary(b.equation_of_state))
+    print(io, "└── equation_of_state: ", summary(b.equation_of_state))
 end
 
 """
@@ -80,7 +80,7 @@ julia> using Oceananigans
 julia> buoyancy = SeawaterBuoyancy(equation_of_state=teos10)
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
-└── equation of state: BoussinesqEquationOfState{Float64}
+└── equation_of_state: BoussinesqEquationOfState{Float64}
 ```
 
 Buoyancy that depends only on salinity with temperature held at 20 degrees Celsius
@@ -90,7 +90,7 @@ julia> salinity_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos10,
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
 ├── constant_temperature: 20
-└── equation of state: BoussinesqEquationOfState{Float64}
+└── equation_of_state: BoussinesqEquationOfState{Float64}
 ```
 
 Buoyancy that depends only on temperature with salinity held at 35 psu
@@ -100,7 +100,7 @@ julia> temperature_dependent_buoyancy = SeawaterBuoyancy(equation_of_state=teos1
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
 ├── constant_salinity: 35
-└── equation of state: BoussinesqEquationOfState{Float64}
+└── equation_of_state: BoussinesqEquationOfState{Float64}
 ```
 """
 function SeawaterBuoyancy(FT = Float64;
diff --git a/src/Models/seawater_density.jl b/src/Models/seawater_density.jl
index f1752ed2c0..67cd252bfc 100644
--- a/src/Models/seawater_density.jl
+++ b/src/Models/seawater_density.jl
@@ -55,16 +55,18 @@ Example
 Compute a density `Field` using the `KernelFunctionOperation` returned from `seawater_density`
 
 ```jldoctest density
-julia> using Oceananigans, Oceananigans.Models, SeawaterPolynomials.TEOS10
+julia> using Oceananigans, SeawaterPolynomials.TEOS10
 
-julia> grid = RectilinearGrid(CPU(), size=(32, 32), x=(0, 2π), y=(0, 2π), topology=(Periodic, Periodic, Flat))
-32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
-├── Periodic x ∈ [3.60072e-17, 6.28319) regularly spaced with Δx=0.19635
-├── Periodic y ∈ [3.60072e-17, 6.28319) regularly spaced with Δy=0.19635
-└── Flat z
+julia> using Oceananigans.Models: seawater_density
 
-julia> tracers = (:S, :T)
-(:S, :T)
+julia> grid = RectilinearGrid(size=100, z=(-1000, 0), topology=(Flat, Flat, Bounded))
+1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
+├── Flat x
+├── Flat y
+└── Bounded  z ∈ [-1000.0, 0.0]   regularly spaced with Δz=10.0
+
+julia> tracers = (:T, :S)
+(:T, :S)
 
 julia> eos = TEOS10EquationOfState()
 BoussinesqEquationOfState{Float64}:
@@ -74,38 +76,44 @@ BoussinesqEquationOfState{Float64}:
 julia> buoyancy = SeawaterBuoyancy(equation_of_state=eos)
 SeawaterBuoyancy{Float64}:
 ├── gravitational_acceleration: 9.80665
-└── equation of state: BoussinesqEquationOfState{Float64}
+└── equation_of_state: BoussinesqEquationOfState{Float64}
 
 julia> model = NonhydrostaticModel(; grid, buoyancy, tracers)
 NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
-├── grid: 32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
+├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
 ├── timestepper: QuasiAdamsBashforth2TimeStepper
-├── tracers: (S, T)
+├── tracers: (T, S)
 ├── closure: Nothing
 ├── buoyancy: SeawaterBuoyancy with g=9.80665 and BoussinesqEquationOfState{Float64} with ĝ = NegativeZDirection()
 └── coriolis: Nothing
 
 julia> set!(model, S = 34.7, T = 0.5)
 
-julia> density_field = Field(seawater_density(model))
-32×32×1 Field{Center, Center, Center} on RectilinearGrid on CPU
-├── grid: 32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
+julia> density_operation = seawater_density(model)
+KernelFunctionOperation at (Center, Center, Center)
+├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
+├── kernel_function: ρ (generic function with 2 methods)
+└── arguments: ("BoussinesqEquationOfState{Float64}", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU", "KernelFunctionOperation at (Center, Center, Center)")
+
+julia> density_field = Field(density_operation)
+1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
+├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
 ├── boundary conditions: FieldBoundaryConditions
-│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
+│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
 ├── operand: KernelFunctionOperation at (Center, Center, Center)
 ├── status: time=0.0
-└── data: 38×38×1 OffsetArray(::Array{Float64, 3}, -2:35, -2:35, 1:1) with eltype Float64 with indices -2:35×-2:35×1:1
+└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
     └── max=0.0, min=0.0, mean=0.0
 
 julia> compute!(density_field)
-32×32×1 Field{Center, Center, Center} on RectilinearGrid on CPU
-├── grid: 32×32×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
+1×1×100 Field{Center, Center, Center} on RectilinearGrid on CPU
+├── grid: 1×1×100 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
 ├── boundary conditions: FieldBoundaryConditions
-│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: ZeroFlux
+│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
 ├── operand: KernelFunctionOperation at (Center, Center, Center)
 ├── status: time=0.0
-└── data: 38×38×1 OffsetArray(::Array{Float64, 3}, -2:35, -2:35, 1:1) with eltype Float64 with indices -2:35×-2:35×1:1
-    └── max=1027.7, min=1027.7, mean=1027.7
+└── data: 1×1×106 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:103) with eltype Float64 with indices 1:1×1:1×-2:103
+    └── max=1032.38, min=1027.73, mean=1030.06
 ```
 
 Values for `temperature`, `salinity` and `geopotential_height` can be passed to