diff --git a/Project.toml b/Project.toml
index 2660e0c..e10171a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,10 +23,12 @@ Variography = "04a0146e-e6df-5636-8d7f-62fa9eb0b20c"
 
 [weakdeps]
 ImageQuilting = "e8712464-036d-575c-85ac-952ae31322ab"
+StratiGraphics = "135379e1-83be-5ae7-9e8e-29dade3dc6c7"
 TuringPatterns = "fde5428d-3bf0-5ade-b94a-d334303c4d77"
 
 [extensions]
 GeoStatsProcessesImageQuiltingExt = "ImageQuilting"
+GeoStatsProcessesStratiGraphicsExt = "StratiGraphics"
 GeoStatsProcessesTuringPatternsExt = "TuringPatterns"
 
 [compat]
@@ -42,6 +44,7 @@ ImageQuilting = "0.22"
 Meshes = "0.35"
 ProgressMeter = "1.9"
 Statistics = "1.9"
+StratiGraphics = "0.7"
 Tables = "1.11"
 TuringPatterns = "0.6"
 Variography = "0.19"
diff --git a/ext/GeoStatsProcessesStratiGraphicsExt.jl b/ext/GeoStatsProcessesStratiGraphicsExt.jl
new file mode 100644
index 0000000..c79ce57
--- /dev/null
+++ b/ext/GeoStatsProcessesStratiGraphicsExt.jl
@@ -0,0 +1,68 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+module GeoStatsProcessesStratiGraphicsExt
+
+using Random
+using Meshes
+
+using StratiGraphics: LandState, Strata, simulate, voxelize
+
+using GeoStatsProcesses: SP, RandSetup
+
+import GeoStatsProcesses: randprep, randsingle
+
+function randprep(::AbstractRNG, process::SP, setup::RandSetup)
+  # retrieve domain info
+  domain = setup.domain
+
+  @assert embeddim(domain) == 3 "solver implemented for 3D domain only"
+
+  pairs = map(setup.varnames) do var
+    # determine initial state
+    state = if isnothing(process.state)
+      nx, ny, _ = size(domain)
+      land = zeros(nx, ny)
+      LandState(land)
+    else
+      process.state
+    end
+
+    # save preprocessed input
+    var => state
+  end
+
+  Dict(pairs)
+end
+
+function randsingle(::AbstractRNG, process::SP, setup::RandSetup, prep)
+  # retrieve domain info
+  domain = setup.domain
+  _, __, nz = size(domain)
+
+  # retrieve process parameters
+  (; environment, stack, nepochs, fillbase, filltop) = process
+
+  pairs = map(setup.varnames) do var
+    # get parameters for the variable
+    state = prep[var]
+
+    # simulate the environment
+    record = simulate(environment, state, nepochs)
+
+    # build stratigraphy
+    strata = Strata(record, stack)
+
+    # return voxel model
+    model = voxelize(strata, nz; fillbase, filltop)
+
+    # flatten result
+    var => vec(model)
+  end
+
+  Dict(pairs)
+end
+
+
+end
diff --git a/src/GeoStatsProcesses.jl b/src/GeoStatsProcesses.jl
index a1b9a0d..3bd4756 100644
--- a/src/GeoStatsProcesses.jl
+++ b/src/GeoStatsProcesses.jl
@@ -33,7 +33,8 @@ include("fft.jl")
 include("lu.jl")
 include("iq.jl")
 include("tp.jl")
+include("sp.jl")
 
-export SPDEGP, SEQ, SGP, FFTGP, LUGP, IQ, TP
+export SPDEGP, SEQ, SGP, FFTGP, LUGP, IQ, TP, SP
 
 end
diff --git a/src/sp.jl b/src/sp.jl
new file mode 100644
index 0000000..f72fbb7
--- /dev/null
+++ b/src/sp.jl
@@ -0,0 +1,12 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+@kwdef struct SP{E,S,ST,N,FB,FT} <: GeoStatsProcess
+  environment::E
+  state::S = nothing
+  stack::ST = :erosional
+  nepochs::N = 10
+  fillbase::FB = NaN
+  filltop::FT = NaN
+end
diff --git a/test/Project.toml b/test/Project.toml
index aefc089..5368a9c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,6 +5,7 @@ ImageQuilting = "e8712464-036d-575c-85ac-952ae31322ab"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StratiGraphics = "135379e1-83be-5ae7-9e8e-29dade3dc6c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TuringPatterns = "fde5428d-3bf0-5ade-b94a-d334303c4d77"
 Variography = "04a0146e-e6df-5636-8d7f-62fa9eb0b20c"
diff --git a/test/runtests.jl b/test/runtests.jl
index e5796cf..62d4220 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,7 @@ using Test
 
 import ImageQuilting
 import TuringPatterns
+using StratiGraphics: SmoothingProcess, Environment, ExponentialDuration
 
 @testset "GeoStatsProcesses.jl" begin
   @testset "FFTGP" begin
@@ -143,5 +144,17 @@ import TuringPatterns
     Random.seed!(2019)
     sdomain = CartesianGrid(200, 200)
     sims = rand(TP(), sdomain, [:z => Float64], 3)
+    @test length(sims) == 3
+    @test size(domain(sims[1])) == (200, 200)
+  end
+
+  @testset "SP" begin
+    rng = MersenneTwister(2019)
+    proc = SmoothingProcess()
+    env = Environment(rng, [proc, proc], [0.5 0.5; 0.5 0.5], ExponentialDuration(rng, 1.0))
+    sdomain = CartesianGrid(50, 50, 20)
+    sims = rand(SP(environment=env), sdomain, [:z => Float64], 3)
+    @test length(sims) == 3
+    @test size(domain(sims[1])) == (50, 50, 20)
   end
 end