diff --git a/src/Denudation/Abstract.jl b/src/Denudation/Abstract.jl
index fb70fcc..90ccf6e 100644
--- a/src/Denudation/Abstract.jl
+++ b/src/Denudation/Abstract.jl
@@ -14,20 +14,9 @@ abstract type DenudationType end
 FIXME Computes the denudation for a single time-step, given denudation parameters `param` and a simulation state `state`. `param` should have a `DenudationType` type and `state` should contain the `height` property and `sealevel`.
 """
 function denudation(input)
-    denudation_mass::Array{typeof(1.0u"m/kyr"),2} = Array{typeof(1.0u"m/kyr")}(undef, input.box.grid_size...)
-
     function (state, water_depth, slope)
-        for idx in CartesianIndices(state.ca)
-            f = state.ca[idx]
-            if f == 0
-                continue
-            end
-
-            if water_depth[idx] >= 0
-                (denudation_mass[idx]) = denudation(input.box, input.denudation, water_depth[idx], slope[idx], input.facies[f])
-            end
-        end
-        return denudation_mass
+        #denu_result = denudation(input.box, input.denudation, water_depth, slope, input.facies)
+        return denudation(input.box, input.denudation, water_depth, slope, input.facies,state)
     end
 end
 
@@ -46,9 +35,13 @@ end
 FIXME
 """
 function redistribution(input)
-    redistribution_mass::Array{typeof(1.0u"m/kyr"),2} = Array{typeof(1.0u"m/kyr")}(undef, input.box.grid_size...)
+        redistribution_mass::Union{Array{typeof(1.0u"m/kyr"),2}, Nothing} = nothing
     function (state, water_depth, denudation_mass)
-        redistribution_mass = redistribution(input.box, input.denudation, denudation_mass, water_depth)
+        redi_result = redistribution(input.box, input.denudation, denudation_mass, water_depth)
+        if redi_result !== nothing
+            redistribution_mass = Array{typeof(1.0u"m/kyr")}(undef, size(state.ca)...)
+            redistribution_mass = redi_result
+        end
         return redistribution_mass
     end
 end
diff --git a/src/Denudation/DissolutionMod.jl b/src/Denudation/DissolutionMod.jl
index 3b91755..a6aa05c 100644
--- a/src/Denudation/DissolutionMod.jl
+++ b/src/Denudation/DissolutionMod.jl
@@ -51,17 +51,27 @@ end
 # ~/~ end
 
 
-function denudation(::Box{BT}, p::Dissolution, water_depth, slope, facies) where {BT<:Boundary}
+function denudation(::Box{BT}, p::Dissolution, water_depth, slope, facies, state) where {BT<:Boundary}
     temp = p.temp ./ u"K"
     precip = p.precip ./u"m"
     pco2 = p.pco2 ./1.0u"atm"
     reactionrate = p.reactionrate ./u"m/yr"
-    return (dissolution(temp, precip, pco2, reactionrate, water_depth, facies))
+    denudation_mass = Array{typeof(1.0u"m/kyr"),2}(undef, size(state.ca)...)
+    for idx in CartesianIndices(state.ca)
+        f = state.ca[idx]
+        if f == 0
+            continue
+        end
+    if water_depth[idx] >= 0
+        denudation_mass[idx] = dissolution(temp, precip, pco2, reactionrate, water_depth[idx], facies[f])
+    end
+    end
+    return denudation_mass
 end
 
 function redistribution(box::Box{BT}, p::Dissolution, denudation_mass, water_depth) where {BT<:Boundary}
-    return denudation_mass .* 0
+    return nothing
 end
 
 end
-# ~/~ end
\ No newline at end of file
+# ~/~ end
diff --git a/src/Denudation/EmpiricalDenudationMod.jl b/src/Denudation/EmpiricalDenudationMod.jl
index b8c3600..4be8630 100644
--- a/src/Denudation/EmpiricalDenudationMod.jl
+++ b/src/Denudation/EmpiricalDenudationMod.jl
@@ -36,13 +36,23 @@ function slope_kernel(w::Any, cellsize::Float64)
 end
 # ~/~ end
 
-function denudation(::Box, p::EmpiricalDenudation, water_depth, slope, facies)
+function denudation(::Box, p::EmpiricalDenudation, water_depth, slope, facies, state)
     precip = p.precip ./ u"m"
-    return empirical_denudation.(precip, slope)
+    denudation_mass = Array{typeof(1.0u"m/kyr"),2}(undef, size(slope)...)
+    for idx in CartesianIndices(state.ca)
+        f = state.ca[idx]
+        if f == 0
+            continue
+        end
+    if water_depth[idx] >= 0
+    denudation_mass[idx] = empirical_denudation.(precip, slope[idx])
+    end
+    end
+    return denudation_mass
 end
 
 function redistribution(::Box, p::EmpiricalDenudation, denudation_mass, water_depth)
-    return denudation_mass .* 0
+    return nothing
 end
 
 end
diff --git a/src/Denudation/NoDenudationMod.jl b/src/Denudation/NoDenudationMod.jl
index ab9ed92..9d84ed5 100644
--- a/src/Denudation/NoDenudationMod.jl
+++ b/src/Denudation/NoDenudationMod.jl
@@ -12,12 +12,12 @@ using Unitful
 
 struct NoDenudation <: DenudationType end
 
-function denudation(box::Box, p::NoDenudation, water_depth::Any, slope, facies)
-    return nothing
+function denudation(box::Box, p::NoDenudation, water_depth::Any, slope, facies, state)
+    return denudation_mass = nothing
 end
 
 function redistribution(box::Box, p::NoDenudation, denudation_mass, water_depth)
-    return denudation_mass .* 0
+    return nothing
 end
 
 end
diff --git a/src/Denudation/PhysicalErosionMod.jl b/src/Denudation/PhysicalErosionMod.jl
index df41bda..9984c8c 100644
--- a/src/Denudation/PhysicalErosionMod.jl
+++ b/src/Denudation/PhysicalErosionMod.jl
@@ -14,7 +14,7 @@ using Unitful
 end
 
 function physical_erosion(slope::Any, inf::Any, erodability::Float64)
-    -1 * -erodability .* (1 - inf) .^ (1 / 3) .* slope .^ (2 / 3)
+    -1 * -erodability .* (1 - inf) .^ (1 / 3) .* slope .^ (2 / 3) .* u"m/kyr"
 end
 
 function redistribution_kernel(w::Array{Float64}, cellsize::Float64)
@@ -63,13 +63,22 @@ function total_mass_redistribution(box::Box{BT}, denudation_mass, water_depth) w
     return mass
 end
 
-function denudation(::Box, p::PhysicalErosion, water_depth::Any, slope, facies)
+function denudation(::Box, p::PhysicalErosion, water_depth::Any, slope, facies,state)
     # This needs transport feature to be merged so that we know the facies type of the
     # top most layer. What follows should still be regarded as pseudo-code.
     # We need to look into this further.
     erodability = p.erodability ./ u"m/yr"
-    denudation_mass = physical_erosion.(slope, facies.infiltration_coefficient, erodability)
-    return (denudation_mass .* u"m/kyr")
+    denudation_mass = Array{typeof(1.0u"m/kyr"),2}(undef, size(slope)...)
+    for idx in CartesianIndices(state.ca)
+        f = state.ca[idx]
+        if f == 0
+            continue
+        end
+    if water_depth[idx] >= 0
+    denudation_mass[idx] = physical_erosion.(slope[idx], facies[f].infiltration_coefficient, erodability)
+    end
+    end
+    return denudation_mass
 end
 
 function redistribution(box::Box{BT}, p::PhysicalErosion, denudation_mass, water_depth) where {BT<:Boundary}
diff --git a/src/Model/WithDenudation.jl b/src/Model/WithDenudation.jl
index 94debba..45e016b 100644
--- a/src/Model/WithDenudation.jl
+++ b/src/Model/WithDenudation.jl
@@ -53,8 +53,8 @@ end
 
 struct OutputFrame
     production::Array{typeof(1.0u"m/Myr"),3}
-    denudation::Array{typeof(1.0u"m/Myr"),2}
-    redistribution::Array{typeof(1.0u"m/Myr"),2}
+    denudation::Union{Array{typeof(1.0u"m/Myr"),2},Nothing}
+    redistribution::Union{Array{typeof(1.0u"m/Myr"),2},Nothing}
 end
 # FIXME: deduplicate
 mutable struct State
@@ -138,8 +138,12 @@ function updater(input)
 
     function update(state, Δ::DenudationFrame)
         # FIXME: implement
+        if Δ.denudation !== nothing
         state.height .+= Δ.denudation .* input.time.Δt
+        end
+        if Δ.redistribution !== nothing
         state.height .-= Δ.redistribution .* input.time.Δt
+        end
         state.time += input.time.Δt
     end
 
@@ -187,24 +191,39 @@ function main(input::Input, output::String)
         attr["subsidence_rate"] = input.subsidence_rate |> in_units_of(u"m/Myr")
         println("Subsidence rate saved successfully.")
 
+        box = input.box
+        results = map(stack_frames, partition(run_model(input, box), input.time.write_interval))
+        testframe = first(results)
 
         n_facies = length(input.facies)
         ds = create_dataset(fid, "sediment", datatype(Float64),
             dataspace(input.box.grid_size..., n_facies, input.time.steps),
             chunk=(input.box.grid_size..., n_facies, 1))
-        denudation = create_dataset(fid, "denudation", datatype(Float64),
-            dataspace(input.box.grid_size..., input.time.steps),
-            chunk=(input.box.grid_size..., 1))
-        redistribution = create_dataset(fid, "redistribution", datatype(Float64),
-            dataspace(input.box.grid_size..., input.time.steps),
-            chunk=(input.box.grid_size..., 1))
-        box = input.box
 
-        results = map(stack_frames, partition(run_model(input, box), input.time.write_interval))
+        denudation = nothing
+        if testframe.denudation !== nothing
+            denudation = create_dataset(fid, "denudation", datatype(Float64),
+                dataspace(input.box.grid_size..., input.time.steps),
+                chunk=(input.box.grid_size..., 1))
+        end
+        
+        redistribution = nothing
+        if testframe.redistribution !== nothing
+            redistribution = create_dataset(fid, "redistribution", datatype(Float64),
+                dataspace(input.box.grid_size..., input.time.steps),
+                chunk=(input.box.grid_size..., 1))
+        end
+
         for (step, frame) in enumerate(take(results, n_writes))
             ds[:, :, :, step] = frame.production |> in_units_of(u"m/Myr")
-            denudation[:, :, step] = frame.denudation |> in_units_of(u"m/kyr")
-            redistribution[:, :, step] = frame.redistribution |> in_units_of(u"m/kyr")
+
+            if denudation !== nothing
+                denudation[:, :, step] = frame.denudation |> in_units_of(u"m/kyr")
+            end
+
+            if redistribution !== nothing
+                redistribution[:, :, step] = frame.redistribution |> in_units_of(u"m/kyr")
+            end
         end
     end
 end