diff --git a/Project.toml b/Project.toml
index b0cf267..02b2a7a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,8 +20,8 @@ Infiltrator = "1.0"
 LineSearches = "7.0"
 NLsolve = "4.5.1"
 PALEOaqchem = "0.3.5"
-PALEOboxes = "0.20.4, 0.21"
-PALEOmodel = "0.15.8"
+PALEOboxes = "0.20.4, 0.21, 0.22"
+PALEOmodel = "0.15.8, 0.16"
 Pkg = "1.0"
 Plots = "1.38"
 SparseArrays = "1.0"
diff --git a/examples/Project.toml b/examples/Project.toml
index 0f11e48..ebbbb0c 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
@@ -16,7 +17,7 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [compat]
 PALEOaqchem = "0.3.10"
-PALEOmodel = "0.15.47"
+PALEOmodel = "0.15.47, 0.16"
 
 [extras]
 PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
diff --git a/src/sediment/SedimentGrid.jl b/src/sediment/SedimentGrid.jl
index bc35a5c..7aedd93 100644
--- a/src/sediment/SedimentGrid.jl
+++ b/src/sediment/SedimentGrid.jl
@@ -235,11 +235,15 @@ function do_sediment_setup_grid(
     end
 
     # attach coordinates to grid for output visualisation etc
-    empty!(rj.domain.grid.z_coords)
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zmid", grid_vars.zmid, PB.get_variable(m, "zmid").attributes))
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zlower", grid_vars.zlower, PB.get_variable(m, "zlower").attributes))
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zupper", grid_vars.zupper, PB.get_variable(m, "zupper").attributes))
-    
+    if isdefined(PB, :set_coordinates!) # PALEOboxes >= 0.22
+        PB.set_coordinates!(rj.domain.grid, "cells", ["zmid", "zlower", "zupper"])
+    else
+        empty!(rj.domain.grid.z_coords)
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zmid", grid_vars.zmid, PB.get_variable(m, "zmid").attributes))
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zlower", grid_vars.zlower, PB.get_variable(m, "zlower").attributes))
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zupper", grid_vars.zupper, PB.get_variable(m, "zupper").attributes))    
+    end
+
     grid_vars.volume .= grid_vars.Abox.*(grid_vars.zupper .- grid_vars.zlower)
     grid_vars.volume_total[] = sum(grid_vars.volume)
 
diff --git a/src/sediment/SedimentTransportDeprecated.jl b/src/sediment/SedimentTransportDeprecated.jl
index 773345e..d27eeea 100644
--- a/src/sediment/SedimentTransportDeprecated.jl
+++ b/src/sediment/SedimentTransportDeprecated.jl
@@ -329,10 +329,14 @@ function do_sediment_setup_grid(
     end
 
     # attach coordinates to grid for output visualisation etc
-    empty!(rj.domain.grid.z_coords)
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zmid", grid_vars.zmid, PB.get_variable(m, "zmid").attributes))
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zlower", grid_vars.zlower, PB.get_variable(m, "zlower").attributes))
-    push!(rj.domain.grid.z_coords, PB.FixedCoord("zupper", grid_vars.zupper, PB.get_variable(m, "zupper").attributes))
+    if isdefined(PB, :set_coordinates!) # PALEOboxes >= 0.22
+        PB.set_coordinates!(rj.domain.grid, "cells", ["zmid", "zlower", "zupper"])
+    else
+        empty!(rj.domain.grid.z_coords)
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zmid", grid_vars.zmid, PB.get_variable(m, "zmid").attributes))
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zlower", grid_vars.zlower, PB.get_variable(m, "zlower").attributes))
+        push!(rj.domain.grid.z_coords, PB.FixedCoord("zupper", grid_vars.zupper, PB.get_variable(m, "zupper").attributes))    
+    end
     
     grid_vars.volume .= grid_vars.Abox.*(grid_vars.zupper .- grid_vars.zlower)
     grid_vars.volume_total[] = sum(grid_vars.volume)