diff --git a/ext/GeoStatsProcessesTuringPatternsExt.jl b/ext/GeoStatsProcessesTuringPatternsExt.jl
index 5542f1b..0dae888 100644
--- a/ext/GeoStatsProcessesTuringPatternsExt.jl
+++ b/ext/GeoStatsProcessesTuringPatternsExt.jl
@@ -21,7 +21,7 @@ function randprep(::AbstractRNG, process::TuringProcess, ::DefaultRandMethod, se
   topo = topology(domain)
 
   # assert grid topology
-  @assert topo isa GridTopology "simulation only defined over grid topology"
+  @assert topo isa GridTopology "process only defined over grid topology"
 
   # retrieve simulation size
   sz = size(topo)
diff --git a/src/field/gaussian.jl b/src/field/gaussian.jl
index 336b15b..e0c1b11 100644
--- a/src/field/gaussian.jl
+++ b/src/field/gaussian.jl
@@ -3,21 +3,18 @@
 # ------------------------------------------------------------------
 
 """
-    GaussianProcess([parameters])
-
-Gaussian process with given variogram and mean.
-
-## Parameters
-
-* `variogram` - Theoretical variogram model (default to `GaussianVariogram()`)
-* `mean`      - Mean field value (default to `0.0`)
+    GaussianProcess(variogram=GaussianVariogram(), mean=0.0)
 
+Gaussian process with given theoretical `variogram` and global `mean`.
 """
-@kwdef struct GaussianProcess{V,T} <: FieldProcess
-  variogram::V = GaussianVariogram()
-  mean::T = 0.0
+struct GaussianProcess{V,T} <: FieldProcess
+  variogram::V
+  mean::T
 end
 
+GaussianProcess(variogram) = GaussianProcess(variogram, 0.0)
+GaussianProcess() = GaussianProcess(GaussianVariogram(), 0.0)
+
 #---------
 # METHODS
 #---------
diff --git a/src/field/gaussian/fft.jl b/src/field/gaussian/fft.jl
index 4c5e38b..62212d4 100644
--- a/src/field/gaussian/fft.jl
+++ b/src/field/gaussian/fft.jl
@@ -3,9 +3,9 @@
 # ------------------------------------------------------------------
 
 """
-    FFTMethod([paramaters])
+    FFTMethod(; [paramaters])
 
-The FFT Gaussian simulation method introduced by Gutjahr 1997.
+The FFT Gaussian process method introduced by Gutjahr 1997.
 The covariance function is perturbed in the frequency domain
 after a fast Fourier transform. White noise is added to the
 phase of the spectrum, and a realization is produced by an
@@ -17,7 +17,7 @@ inverse Fourier transform.
 * `neighborhood` - Search neighborhood (default to `nothing`)
 * `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)
 
-### References
+## References
 
 * Gutjahr 1997. [General joint conditional simulations using a fast
   Fourier transform method](https://link.springer.com/article/10.1007/BF02769641)
diff --git a/src/field/gaussian/lu.jl b/src/field/gaussian/lu.jl
index ca2c48d..4f14a60 100644
--- a/src/field/gaussian/lu.jl
+++ b/src/field/gaussian/lu.jl
@@ -3,11 +3,11 @@
 # ------------------------------------------------------------------
 
 """
-    LUMethod([paramaters])
+    LUMethod(; [paramaters])
 
-The LU Gaussian simulation method introduced by Alabert 1987.
+The LU Gaussian process method introduced by Alabert 1987.
 The full covariance matrix is built to include all locations
-of the simulation domain, and samples from the multivariate
+of the process domain, and samples from the multivariate
 Gaussian are drawn via LU factorization.
 
 ## Parameters
@@ -16,7 +16,7 @@ Gaussian are drawn via LU factorization.
 * `correlation`   - Correlation coefficient between two covariates (default to `0`)
 * `init`          - Data initialization method (default to `NearestInit()`)
 
-### References
+## References
 
 * Alabert 1987. [The practice of fast conditional simulations
   through the LU decomposition of the covariance matrix]
diff --git a/src/field/gaussian/seq.jl b/src/field/gaussian/seq.jl
index e815ca3..cc633cd 100644
--- a/src/field/gaussian/seq.jl
+++ b/src/field/gaussian/seq.jl
@@ -3,9 +3,9 @@
 # ------------------------------------------------------------------
 
 """
-    SEQMethod([paramaters])
+    SEQMethod(; [paramaters])
 
-The sequential simulation method introduced by Gomez-Hernandez 1993.
+The sequential process method introduced by Gomez-Hernandez 1993.
 It traverses all locations of the geospatial domain according to a path,
 approximates the conditional Gaussian distribution within a neighborhood
 using Kriging, and assigns a value to the center of the neighborhood by
@@ -13,25 +13,25 @@ sampling from this distribution.
 
 ## Parameters
 
-* `path`         - Simulation path (default to `LinearPath()`)
+* `path`         - Process path (default to `LinearPath()`)
 * `minneighbors` - Minimum number of neighbors (default to `1`)
 * `maxneighbors` - Maximum number of neighbors (default to `10`)
 * `neighborhood` - Search neighborhood (default to `nothing`)
 * `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)
 * `init`         - Data initialization method (default to `NearestInit()`)
 
-For each location in the simulation `path`, a maximum number of
+For each location in the process `path`, a maximum number of
 neighbors `maxneighbors` is used to fit the conditional Gaussian
 distribution. The neighbors are searched according to a `neighborhood`.
 
-### References
+## References
 
 * Gomez-Hernandez & Journel 1993. [Joint Sequential Simulation of
   MultiGaussian Fields](https://link.springer.com/chapter/10.1007/978-94-011-1739-5_8)
 
 ### Notes
 
-* This method is very sensitive to the simulation path and number of
+* This method is very sensitive to the process path and number of
   samples. Care must be taken to make sure that neighborhoods have
   enough samples for the geostatistical model (e.g. Kriging).
 """
@@ -49,7 +49,7 @@ function randprep(::AbstractRNG, process::GaussianProcess, method::SEQMethod, se
   (; path, minneighbors, maxneighbors, neighborhood, distance, init) = method
   probmodel = GeoStatsModels.SimpleKriging(variogram, mean)
   marginal = Normal(mean, √sill(variogram))
-  SequentialProcess(; probmodel, marginal, path, minneighbors, maxneighbors, neighborhood, distance, init)
+  SequentialProcess(probmodel, marginal; path, minneighbors, maxneighbors, neighborhood, distance, init)
 end
 
 function randsingle(rng::AbstractRNG, ::GaussianProcess, ::SEQMethod, setup::RandSetup, prep)
diff --git a/src/field/lindgren.jl b/src/field/lindgren.jl
index c67e575..1a4678c 100644
--- a/src/field/lindgren.jl
+++ b/src/field/lindgren.jl
@@ -3,30 +3,30 @@
 # ------------------------------------------------------------------
 
 """
-    LindgrenProcess([paramaters])
+    LindgrenProcess(range=1.0, sill=1.0)
 
-The SPDE process introduced by Lindgren 2011. It relies on a
-discretization of the Laplace-Beltrami operator on meshes and
-is adequate for highly curved domains (e.g. surfaces).
+Lindgren process with given `range` (correlation length)
+and `sill` (total variance) as described in Lindgren 2011.
 
-## Parameters
+The process relies relies on a discretization of the Laplace-Beltrami 
+operator on meshes and is adequate for highly curved domains (e.g. surfaces).
 
-* `range` - Range or correlation length (default to `1.0`)
-* `sill`  - Sill or total variance (default to `1.0`)
-
-### References
+## References
 
 * Lindgren et al. 2011. [An explicit link between Gaussian fields and
   Gaussian Markov random fields: the stochastic partial differential
   equation approach](https://rss.onlinelibrary.wiley.com/doi/10.1111/j.1467-9868.2011.00777.x)
 """
-@kwdef struct LindgrenProcess <: FieldProcess
-  range::Float64 = 1.0
-  sill::Float64 = 1.0
+struct LindgrenProcess <: FieldProcess
+  range::Float64
+  sill::Float64
 end
 
+LindgrenProcess(range) = LindgrenProcess(range, 1.0)
+LindgrenProcess() = LindgrenProcess(1.0, 1.0)
+
 function randprep(::AbstractRNG, process::LindgrenProcess, ::DefaultRandMethod, setup::RandSetup)
-  isnothing(setup.geotable) || @error "conditional simulation is not implemented"
+  isnothing(setup.geotable) || @error "conditional process is not implemented"
 
   # retrieve sill and range
   𝓁 = process.range
diff --git a/src/field/quilting.jl b/src/field/quilting.jl
index 49c8c0e..444e797 100644
--- a/src/field/quilting.jl
+++ b/src/field/quilting.jl
@@ -3,21 +3,15 @@
 # ------------------------------------------------------------------
 
 """
-    QuiltingProcess([paramaters])
+    QuiltingProcess(trainimg, tilesize; [paramaters])
 
-Image quilting process as described in Hoffimann et al. 2017.
+Image quilting process with training image `trainimg` and tile size `tilesize`
+as described in Hoffimann et al. 2017.
 
 ## Parameters
 
-### Required
-
-* `trainimg` - Training image from which to extract tiles
-* `tilesize` - Tuple with tile size for each dimension
-
-### Optional
-
 * `overlap`  - Overlap size (default to (1/6, 1/6, ..., 1/6))
-* `path`     - Simulation path (`:raster` (default), `:dilation`, or `:random`)
+* `path`     - Process path (`:raster` (default), `:dilation`, or `:random`)
 * `inactive` - Vector of inactive voxels (i.e. `CartesianIndex`) in the grid
 * `soft`     - A pair `(data,dataTI)` of geospatial data objects (default to `nothing`)
 * `tol`      - Initial relaxation tolerance in (0,1] (default to `0.1`)
@@ -25,16 +19,30 @@ Image quilting process as described in Hoffimann et al. 2017.
 
 ## References
 
-* Hoffimann et al 2017. *Stochastic simulation by image quilting of process-based geological models.*
-* Hoffimann et al 2015. *Geostatistical modeling of evolving landscapes by means of image quilting.*
+* Hoffimann et al 2017. [Stochastic simulation by image quilting
+  of process-based geological models](https://www.sciencedirect.com/science/article/abs/pii/S0098300417301139)
+
+* Hoffimann et al 2015. [Geostatistical modeling of evolving landscapes
+  by means of image quilting](https://www.researchgate.net/publication/295902985_Geostatistical_Modeling_of_Evolving_Landscapes_by_Means_of_Image_Quilting)
 """
-@kwdef struct QuiltingProcess{TR,TS,O,P,IN,S,T,I} <: FieldProcess
+struct QuiltingProcess{TR,TS,O,P,IN,S,T,I} <: FieldProcess
   trainimg::TR
   tilesize::TS
-  overlap::O = nothing
-  path::P = :raster
-  inactive::IN = nothing
-  soft::S = nothing
-  tol::T = 0.1
-  init::I = NearestInit()
+  overlap::O
+  path::P
+  inactive::IN
+  soft::S
+  tol::T
+  init::I
 end
+
+QuiltingProcess(
+  trainimg,
+  tilesize;
+  overlap=nothing,
+  path=:raster,
+  inactive=nothing,
+  soft=nothing,
+  tol=0.1,
+  init=NearestInit()
+) = QuiltingProcess(trainimg, tilesize, overlap, path, inactive, soft, tol, init)
diff --git a/src/field/sequential.jl b/src/field/sequential.jl
index 2fde107..218710c 100644
--- a/src/field/sequential.jl
+++ b/src/field/sequential.jl
@@ -2,39 +2,48 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
-
 """
-    SequentialProcess([paramaters])
+    SequentialProcess(probmodel, marginal; [paramaters])
 
-A sequential simulation process.
+A sequential process with given conditional distribution model `probmodel`
+and `marginal` distribution.
 
-For each location in the simulation `path`, a maximum number
+For each location in the process `path`, a maximum number
 of neighbors `maxneighbors` is used to fit a distribution.
 The neighbors are searched according to a `neighborhood`,
 and in case there are none, use a `marginal` distribution.
 
 ## Parameters
 
-* `probmodel`    - Conditional distribution model from GeoStatsModels.jl
-* `marginal`     - Marginal distribution from Distributions.jl
-* `path`         - Simulation path (default to `LinearPath()`)
+* `path`         - Process path (default to `LinearPath()`)
 * `minneighbors` - Minimum number of neighbors (default to `1`)
 * `maxneighbors` - Maximum number of neighbors (default to `10`)
 * `neighborhood` - Search neighborhood (default to `nothing`)
 * `distance`     - Distance used to find nearest neighbors (default to `Euclidean()`)
 * `init`         - Data initialization method (default to `NearestInit()`)
 """
-@kwdef struct SequentialProcess{M,MD,P,N,D,I} <: FieldProcess
+struct SequentialProcess{M,MD,P,N,D,I} <: FieldProcess
   probmodel::M
   marginal::MD
-  path::P = LinearPath()
-  minneighbors::Int = 1
-  maxneighbors::Int = 10
-  neighborhood::N = nothing
-  distance::D = Euclidean()
-  init::I = NearestInit()
+  path::P
+  minneighbors::Int
+  maxneighbors::Int
+  neighborhood::N
+  distance::D
+  init::I
 end
 
+SequentialProcess(
+  probmodel,
+  marginal;
+  path=LinearPath(),
+  minneighbors=1,
+  maxneighbors=10,
+  neighborhood=nothing,
+  distance=Euclidean(),
+  init=NearestInit()
+) = SequentialProcess(probmodel, marginal, path, minneighbors, maxneighbors, neighborhood, distance, init)
+
 function randprep(::AbstractRNG, process::SequentialProcess, ::DefaultRandMethod, setup::RandSetup)
   # retrieve domain info
   domain = setup.domain
diff --git a/src/field/strata.jl b/src/field/strata.jl
index 5e9b632..c1ac694 100644
--- a/src/field/strata.jl
+++ b/src/field/strata.jl
@@ -3,29 +3,32 @@
 # ------------------------------------------------------------------
 
 """
-    StrataProcess([paramaters])
+    StrataProcess(environment; [paramaters])
 
-A process to simulate 3D stratigraphic models.
+Strata process with given geological `environment`
+as described in Hoffimann 2018.
 
 ## Parameters
 
-* `environment` - Geological environment
 * `state`       - Initial geological state
 * `stack`       - Stacking scheme (:erosional or :depositional)
 * `nepochs`     - Number of epochs (default to 10)
 * `fillbase`    - Fill value for the bottom layer (default to `NaN`)
 * `filltop`     - Fill value for the top layer (default to `NaN`)
 
-### References
+## References
 
-Hoffimann 2018. *Morphodynamic analysis and statistical
-synthesis of geormorphic data.*
+* Hoffimann 2018. [Morphodynamic analysis and statistical synthesis of
+  geormorphic data](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JF005245)
 """
-@kwdef struct StrataProcess{E,S,ST,N,FB,FT} <: FieldProcess
+struct StrataProcess{E,S,ST,N,FB,FT} <: FieldProcess
   environment::E
-  state::S = nothing
-  stack::ST = :erosional
-  nepochs::N = 10
-  fillbase::FB = NaN
-  filltop::FT = NaN
+  state::S
+  stack::ST
+  nepochs::N
+  fillbase::FB
+  filltop::FT
 end
+
+StrataProcess(environment; state=nothing, stack=:erosional, nepochs=10, fillbase=NaN, filltop=NaN) =
+  StrataProcess(environment, state, stack, nepochs, fillbase, filltop)
diff --git a/src/field/turing.jl b/src/field/turing.jl
index 23ff93c..9a7aca1 100644
--- a/src/field/turing.jl
+++ b/src/field/turing.jl
@@ -3,7 +3,7 @@
 # ------------------------------------------------------------------
 
 """
-    TuringProcess([paramaters])
+    TuringProcess(; [paramaters])
 
 Turing process as described in Turing 1952.
 
@@ -14,13 +14,15 @@ Turing process as described in Turing 1952.
 * `edge`   - edge condition (default to `nothing`)
 * `iter`   - number of iterations (default to `100`)
 
-### References
+## References
 
-Turing 1952. *The chemical basis of morphogenesis.*
+* Turing 1952. [The chemical basis of morphogenesis](https://royalsocietypublishing.org/doi/10.1098/rstb.1952.0012)
 """
-@kwdef struct TuringProcess{P,B,E,I} <: FieldProcess
-  params::P = nothing
-  blur::B = nothing
-  edge::E = nothing
-  iter::I = 100
+struct TuringProcess{P,B,E,I} <: FieldProcess
+  params::P
+  blur::B
+  edge::E
+  iter::I
 end
+
+TuringProcess(; params=nothing, blur=nothing, edge=nothing, iter=100) = TuringProcess(params, blur, edge, iter)
diff --git a/test/processes.jl b/test/processes.jl
index 03c2624..a9ee2bc 100644
--- a/test/processes.jl
+++ b/test/processes.jl
@@ -22,14 +22,14 @@
         # isotropic simulation
         Random.seed!(2019)
         dom = CartesianGrid(100, 100)
-        process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
+        process = GaussianProcess(GaussianVariogram(range=10.0))
         method = FFTMethod()
         sims = rand(process, dom, [:z => Float64], 3, method)
 
         # anisotropic simulation
         Random.seed!(2019)
         dom = CartesianGrid(100, 100)
-        process = GaussianProcess(variogram=GaussianVariogram(MetricBall((20.0, 5.0))))
+        process = GaussianProcess(GaussianVariogram(MetricBall((20.0, 5.0))))
         method = FFTMethod()
         sims = rand(process, dom, [:z => Float64], 3, method)
 
@@ -37,7 +37,7 @@
         Random.seed!(2022)
         grid = CartesianGrid(100, 100)
         vgrid = view(grid, 1:5000)
-        process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
+        process = GaussianProcess(GaussianVariogram(range=10.0))
         method = FFTMethod()
         sim = rand(process, vgrid, [:z => Float64], method)
         @test domain(sim) == vgrid
@@ -49,7 +49,7 @@
         coords = [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)]
         samples = georef(table, coords)
         sdomain = CartesianGrid(100, 100)
-        process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
+        process = GaussianProcess(GaussianVariogram(range=10.0))
         method = FFTMethod(maxneighbors=3)
         sim = rand(process, sdomain, samples, method)
       end
@@ -59,7 +59,7 @@
         𝒟 = CartesianGrid((100, 100), (0.5, 0.5), (1.0, 1.0))
         N = 3
 
-        process = GaussianProcess(variogram=SphericalVariogram(range=35.0))
+        process = GaussianProcess(SphericalVariogram(range=35.0))
         method = SEQMethod(neighborhood=MetricBall(30.0), maxneighbors=3)
 
         Random.seed!(2017)
@@ -82,7 +82,7 @@
         # conditional simulation
         # ----------------------
         rng = MersenneTwister(123)
-        process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
+        process = GaussianProcess(SphericalVariogram(range=10.0))
         method = LUMethod()
         sims = rand(rng, process, 𝒟, 𝒮, 2, method)
 
@@ -90,7 +90,7 @@
         # unconditional simulation
         # ------------------------
         rng = MersenneTwister(123)
-        process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
+        process = GaussianProcess(SphericalVariogram(range=10.0))
         method = LUMethod()
         sims = rand(rng, process, 𝒟, [:z => Float64], 2, method)
 
@@ -99,7 +99,7 @@
         # -------------
         𝒟 = CartesianGrid(500)
         rng = MersenneTwister(123)
-        process = GaussianProcess(variogram=(SphericalVariogram(range=10.0), GaussianVariogram(range=10.0)))
+        process = GaussianProcess((SphericalVariogram(range=10.0), GaussianVariogram(range=10.0)))
         method = LUMethod(correlation=0.95)
         sim = rand(rng, process, 𝒟, [:a => Float64, :b => Float64], method)
 
@@ -108,7 +108,7 @@
         # -----------
         𝒟 = CartesianGrid(100, 100)
         rng = MersenneTwister(123)
-        process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
+        process = GaussianProcess(GaussianVariogram(range=10.0))
         method = LUMethod()
         sims = rand(rng, process, 𝒟, [:z => Float64], 3, method)
 
@@ -118,7 +118,7 @@
         𝒟 = CartesianGrid(100, 100)
         rng = MersenneTwister(123)
         ball = MetricBall((20.0, 5.0))
-        process = GaussianProcess(variogram=GaussianVariogram(ball))
+        process = GaussianProcess(GaussianVariogram(ball))
         method = LUMethod()
         sims = rand(rng, process, 𝒟, [:z => Float64], 3, method)
 
@@ -127,7 +127,7 @@
         # ---------------------
         𝒟 = CartesianGrid(100)
         rng = MersenneTwister(123)
-        process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
+        process = GaussianProcess(SphericalVariogram(range=10.0))
         method1 = LUMethod(factorization=lu)
         method2 = LUMethod(factorization=cholesky)
         sim1 = rand(rng, process, 𝒟, 𝒮, 2, method1)
@@ -135,11 +135,11 @@
 
         # throws
         𝒟 = CartesianGrid(100, 100)
-        process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
+        process = GaussianProcess(GaussianVariogram(range=10.0))
         method = LUMethod()
         # only 1 or 2 variables can be simulated simultaneously
         @test_throws AssertionError rand(process, 𝒟, [:a => Float64, :b => Float64, :c => Float64], method)
-        process = GaussianProcess(variogram=(GaussianVariogram(range=10.0),))
+        process = GaussianProcess((GaussianVariogram(range=10.0),))
         # the number of parameters must be equal to the number of variables
         @test_throws AssertionError rand(process, 𝒟, [:a => Float64, :b => Float64], method)
       end
@@ -152,7 +152,7 @@
       rng = MersenneTwister(2017)
       trainimg = geostatsimage("Strebelle")
       inactive = [CartesianIndex(i, j) for i in 1:30 for j in 1:30]
-      process = QuiltingProcess(trainimg=trainimg, tilesize=(30, 30), inactive=inactive)
+      process = QuiltingProcess(trainimg, (30, 30); inactive)
 
       sims = rand(rng, process, sdomain, sdata, 3)
       @test length(sims) == 3
@@ -172,7 +172,7 @@
       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(StrataProcess(environment=env), sdomain, [:z => Float64], 3)
+      sims = rand(StrataProcess(env), sdomain, [:z => Float64], 3)
       @test length(sims) == 3
       @test size(domain(sims[1])) == (50, 50, 20)
     end