diff --git a/Project.toml b/Project.toml
index e2fb067..c89c051 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,7 +21,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 CUDA = "1.3, 2, 3.7.1, 4, 5"
 DocStringExtensions = "0.8, 0.9"
 FourierFlows = "0.10.5"
-GeophysicalFlows = "0.15.4"
+GeophysicalFlows = "0.16"
 JLD2 = "0.1, 0.2, 0.3, 0.4"
 LinearAlgebra = "1.6"
 Random = "1.6"
diff --git a/examples/turbulent_advection-diffusion.jl b/examples/turbulent_advection-diffusion.jl
index e7081b8..b9d0cad 100644
--- a/examples/turbulent_advection-diffusion.jl
+++ b/examples/turbulent_advection-diffusion.jl
@@ -42,18 +42,18 @@ L = 2π                   # domain size
 β = 5                    # the y-gradient of planetary PV
 
 nlayers = 2              # number of layers
-f₀, g = 1, 1             # Coriolis parameter and gravitational constant
+f₀ = 1                   # Coriolis parameter
 H = [0.2, 0.8]           # the rest depths of each layer
-ρ = [4.0, 5.0]           # the density of each layer
+b = [-1.0, -1.2]         # Boussinesq buoyancy of each layer
 
- U = zeros(nlayers)      # the imposed mean zonal flow in each layer
- U[1] = 1.0
- U[2] = 0.0
+U = zeros(nlayers)       # the imposed mean zonal flow in each layer
+U[1] = 1.0
+U[2] = 0.0
 nothing # hide
 
 # ### `MultiLayerQG.Problem` setup, shortcuts and initial conditions
 MQGprob = MultiLayerQG.Problem(nlayers, dev;
-                               nx=n, Lx=L, f₀, g, H, ρ, U, μ, β,
+                               nx=n, Lx=L, f₀, H, b, U, μ, β,
                                dt, stepper, aliased_fraction=0)
 
 nx, ny = MQGprob.grid.nx, MQGprob.grid.ny
diff --git a/test/test_traceradvectiondiffusion.jl b/test/test_traceradvectiondiffusion.jl
index 421a5f7..91fc052 100644
--- a/test/test_traceradvectiondiffusion.jl
+++ b/test/test_traceradvectiondiffusion.jl
@@ -305,20 +305,20 @@ function test_diffusion_multilayerqg(stepper, dt, tfinal, dev::Device=CPU())
   nx = 128
   Lx = 2π
 
-  μ = 0                 
-  β = 0                    
+  μ = 0
+  β = 0
 
-  nlayers = 2              
-  f₀, g = 1, 1             
+  nlayers = 2
+  f₀ = 1
   H = [0.2, 0.8]          
-  ρ = [4.0, 5.0]          
+  b = [-1.0, -1.2]
 
   U = zeros(nlayers) 
   U[1] = 0.0
   U[2] = 0.0
 
   MQGprob = MultiLayerQG.Problem(nlayers, dev;
-                                 nx, Lx, f₀, g, H, ρ, U, μ, β, dt,
+                                 nx, Lx, f₀, H, b, U, μ, β, dt,
                                  stepper="FilteredRK4", aliased_fraction=0)
   grid = MQGprob.grid
   q₀ = zeros(dev, eltype(grid), (grid.nx, grid.ny, nlayers))