diff --git a/docs/source/reference/model/model_hooks.rst b/docs/source/reference/model/model_hooks.rst
index 99584697..348ab218 100644
--- a/docs/source/reference/model/model_hooks.rst
+++ b/docs/source/reference/model/model_hooks.rst
@@ -57,4 +57,6 @@ A complete list of behavior-modifying arrays in the model follows:
     :widths: 20, 50, 10
 
     `mod_water_weight`, modifies the neighbor-weighting of water parcels during routing according to ``(depth * mod_water_weight)**theta_water``, 1
+    `mod_sed_weight`, modifies the neighbor-weighting of the sediment parcel routing according to ``(depth * mod_sed_weight)**theta_sed``, 1
+    `mod_erosion`, linearly modifies the "erodibility" of cells according to ``mod_erosion * Vp_sed * (U_loc**beta - U_ero**beta) / U_ero**beta``, 1
 
diff --git a/pyDeltaRCM/init_tools.py b/pyDeltaRCM/init_tools.py
index 54867444..a4071000 100644
--- a/pyDeltaRCM/init_tools.py
+++ b/pyDeltaRCM/init_tools.py
@@ -257,7 +257,7 @@ def create_other_variables(self) -> None:
 
         Creates variables for model implementation, from specified boundary
         condition variables. This method is run during initial model
-        instantition. Internally, this method calls :obj:`set_constants` and
+        instantiation. Internally, this method calls :obj:`set_constants` and
         :obj:`create_boundary_conditions`.
 
         .. note::
@@ -452,7 +452,7 @@ def create_domain(self) -> None:
 
             If you need to modify the model domain after it has been created,
             it is probably safe to modify attributes directly, but take care
-            to ensure any dependent fields are also approrpriately changed.
+            to ensure any dependent fields are also appropriately changed.
         """
         _msg = 'Creating model domain'
         self.log_info(_msg, verbosity=1)
@@ -499,6 +499,8 @@ def create_domain(self) -> None:
         
         # arrays acting as modifying hooks
         self.mod_water_weight = np.ones_like(self.depth)
+        self.mod_sed_weight = np.ones_like(self.depth)
+        self.mod_erosion = np.ones_like(self.depth)
 
         # ---- domain ----
         cell_land = -2
@@ -571,7 +573,7 @@ def init_sediment_routers(self) -> None:
                                        self.iwalk_flat, self.jwalk_flat,
                                        self.distances_flat, self.dry_depth,
                                        self._lambda, self._beta,  self.stepmax,
-                                       self.theta_mud)
+                                       self.theta_mud, self.mod_erosion)
         # initialize the SandRouter object
         self._sr = sed_tools.SandRouter(self._dt, self._dx, self.Vp_sed,
                                         self.u_max, self.qs0, self._u0,
@@ -581,7 +583,7 @@ def init_sediment_routers(self) -> None:
                                         self.iwalk_flat, self.jwalk_flat,
                                         self.distances_flat, self.dry_depth,
                                         self._beta, self.stepmax,
-                                        self.theta_sand)
+                                        self.theta_sand, self.mod_erosion)
 
     def init_output_file(self) -> None:
         """Creates a netCDF file to store output grids.
@@ -779,7 +781,7 @@ def load_checkpoint(self, defer_output: bool = False) -> None:
 
         As a standard user, you should not need to worry about any of these
         pathways or options. However, if you are developing pyDeltaRCM or
-        customizing the model in any way that involves loadind from
+        customizing the model in any way that involves loading from
         checkpoints, you should be aware of these pathways.
 
         For example, loading from checkpoint will succeed if no netCDF4 file
@@ -789,7 +791,7 @@ def load_checkpoint(self, defer_output: bool = False) -> None:
 
         .. important::
 
-            If you are customing the model and intend to use checkpointing and
+            If you are customizing the model and intend to use checkpointing and
             the :obj:`Preprocessor` parallel infrastructure, be sure that
             parameter :obj:`defer_output` is `True` until the
             :obj:`load_checkpoint` method can be called from the thread the
diff --git a/pyDeltaRCM/iteration_tools.py b/pyDeltaRCM/iteration_tools.py
index cda5d0c5..3f5fb950 100644
--- a/pyDeltaRCM/iteration_tools.py
+++ b/pyDeltaRCM/iteration_tools.py
@@ -235,8 +235,8 @@ def compute_sand_frac(self) -> None:
     def save_grids_and_figs(self) -> None:
         """Save grids and figures.
 
-        Save grids and/or plots of specified variables (``eta``, `discharge``,
-        ``velocity``, ``depth``, and ``stage``, depending on configuration of
+        Save grids and/or plots of specified variables (``eta``, ``discharge``,
+        ``velocity``, ``depth``, and ``stage``), depending on configuration of
         the relevant flags in the YAML configuration file.
 
         .. note:
diff --git a/pyDeltaRCM/model.py b/pyDeltaRCM/model.py
index 08eeb7c8..4aa5f68d 100644
--- a/pyDeltaRCM/model.py
+++ b/pyDeltaRCM/model.py
@@ -48,7 +48,7 @@ def __init__(self, input_file=None, defer_output: bool = False, **kwargs: Any) -
             Whether to create the output netCDF file during initialization. In
             most cases, this can be ignored and left to default value `False`.
             However, for parallel simulations it may be necessary to defer the
-            NetCDF file creation unitl the simualtion is assigned to the core
+            NetCDF file creation until the simulation is assigned to the core
             it will compute on, so that the DeltaModel object remains
             pickle-able. Note, you will need to manually trigger the
             :obj:`init_output_file`, :obj:`output_data`, and
@@ -59,7 +59,7 @@ def __init__(self, input_file=None, defer_output: bool = False, **kwargs: Any) -
             be passed as a keyword argument to the instantiation of the
             DeltaModel. Keywords will override the specification of any value
             in the YAML file. This functionality is intended for advanced use
-            cases, and should not be preffered to specifying all inputs in a
+            cases, and should not be preferred to specifying all inputs in a
             YAML configuration file.
 
         Returns
diff --git a/pyDeltaRCM/preprocessor.py b/pyDeltaRCM/preprocessor.py
index aef60fb5..e00e61ff 100644
--- a/pyDeltaRCM/preprocessor.py
+++ b/pyDeltaRCM/preprocessor.py
@@ -1127,7 +1127,7 @@ def preprocessor_wrapper() -> None:
 def _write_yaml_config_to_file(_config, _path) -> None:
     """Write a config to file in output folder.
 
-    Write the entire yaml configuation for the configured job out to a
+    Write the entire yaml configuration for the configured job out to a
     file in the job output foler.
 
     .. note::
diff --git a/pyDeltaRCM/sed_tools.py b/pyDeltaRCM/sed_tools.py
index 9a0212c7..17ccd157 100644
--- a/pyDeltaRCM/sed_tools.py
+++ b/pyDeltaRCM/sed_tools.py
@@ -129,7 +129,7 @@ def route_all_sand_parcels(self) -> None:
         self._sr.run(start_indices, self.eta, self.stage, self.depth,
                      self.cell_type, self.uw, self.ux, self.uy,
                      self.Vp_dep_mud, self.Vp_dep_sand,
-                     self.qw, self.qx, self.qy, self.qs)
+                     self.qw, self.qx, self.qy, self.qs, self.mod_sed_weight)
 
         # These are the variables updated at the end of the `SandRouter`. If
         # you attempt to drop in a replacement SandRouter, you will need to
@@ -184,7 +184,7 @@ def route_all_mud_parcels(self) -> None:
         self._mr.run(start_indices, self.eta, self.stage, self.depth,
                      self.cell_type, self.uw, self.ux, self.uy,
                      self.Vp_dep_mud, self.Vp_dep_sand,
-                     self.qw, self.qx, self.qy)
+                     self.qw, self.qx, self.qy, self.mod_sed_weight)
 
         # These are the variables updated at the end of the `MudRouter`. If
         # you attempt to drop in a replacement MudRouter, you will need to
@@ -225,7 +225,8 @@ def topo_diffusion(self) -> None:
 
 @njit
 def _get_weight_at_cell_sediment(ind, weight_int, depth_nbrs, ct_nbrs,
-                                 dry_depth: float, theta: float, distances_flat: int):
+                                 dry_depth: float, theta: float,
+                                 distances_flat: int, mod_sed_weight):
     """Get neighbor weight array for sediment routing.
 
     .. todo::
@@ -242,7 +243,7 @@ def _get_weight_at_cell_sediment(ind, weight_int, depth_nbrs, ct_nbrs,
     weight_int[ctr] = 0
 
     weight = np.copy(weight_int)  # no gamma weighting here
-    weight = (depth_nbrs ** theta) * weight
+    weight = ((depth_nbrs * mod_sed_weight) ** theta) * weight
 
     # ALWAYS disallow the choice to not move
     weight[ctr] = 0
@@ -298,7 +299,9 @@ def _get_weight_at_cell_sediment(ind, weight_int, depth_nbrs, ct_nbrs,
           ('Vp_res', float32), ('Vp_dep_mud', float32[:, :]),
           ('Vp_dep_sand', float32[:, :]),
           ('U_dep_mud', float32), ('U_ero_mud', float32),
-          ('U_ero_sand', float32)]
+          ('U_ero_sand', float32), ('mod_erosion', float32[:, :]),
+          ('mod_sed_weight', float32[:, :]),
+          ('pad_mod_sed_weight', float32[:, :])]
 
 
 class BaseRouter:
@@ -336,6 +339,8 @@ def _choose_next_location(self, px: int, py: int) -> Tuple[int, int, float]:
             px:px + 3, py:py + 3]
         cell_type_ind = self.pad_cell_type[
             px:px + 3, py:py + 3]
+        sed_weight_nbrs = self.pad_mod_sed_weight[
+            px:px + 3, py:py + 3]
 
         _, weight_int = shared_tools.get_weight_sfc_int(
             self.stage[px, py], stage_nbrs.ravel(), self.qx[px, py],
@@ -350,7 +355,8 @@ def _choose_next_location(self, px: int, py: int) -> Tuple[int, int, float]:
 
         weights = _get_weight_at_cell_sediment(
             (px, py), weight_int, depth_nbrs.ravel(),
-            cell_type_ind.ravel(), self.dry_depth, self.theta_sed, self.distances_flat)
+            cell_type_ind.ravel(), self.dry_depth, self.theta_sed,
+            self.distances_flat, sed_weight_nbrs.ravel())
 
         new_cell = shared_tools.random_pick(weights)
 
@@ -435,7 +441,9 @@ def _update_fields(self, Vp_change: float, px: int, py: int) -> None:
             self.ux[px, py] = 0
             self.uy[px, py] = 0
 
-    def _compute_Vp_ero(self, Vp_sed: float, U_loc: float, U_ero: float, beta: float) -> float:
+    def _compute_Vp_ero(self, Vp_sed: float, U_loc: float,
+                        U_ero: float, beta: float,
+                        mod_erosion: float) -> float:
         """Compute volume erorded based on velocity.
 
         The volume of sediment eroded depends on the local flow velocity
@@ -459,12 +467,15 @@ def _compute_Vp_ero(self, Vp_sed: float, U_loc: float, U_ero: float, beta: float
         beta : :obj:`float`
             unknown.
 
+        mod_erosion : :obj:`float`
+            Local linear modifier for erosion.
+
         Returns
         -------
         Vp_ero : :obj:`float`
             Volume of eroded sediment.
         """
-        return (Vp_sed * (U_loc**beta - U_ero**beta) /
+        return (mod_erosion * Vp_sed * (U_loc**beta - U_ero**beta) /
                 U_ero**beta)
 
     def _limit_Vp_change(self, Vp, stage, eta, dx: float, dep_ero: int):
@@ -507,7 +518,7 @@ class SandRouter(BaseRouter):
     def __init__(self, _dt: float, dx: float, Vp_sed, u_max: float, qs0: float,
                  u0, U_ero_sand, f_bedload, ivec_flat, jvec_flat, iwalk_flat,
                  jwalk_flat, distances_flat, dry_depth: float, beta: float,
-                 stepmax, theta_sed: float) -> None:
+                 stepmax, theta_sed: float, mod_erosion) -> None:
 
         self._dt = _dt
         self._dx = dx
@@ -527,10 +538,11 @@ def __init__(self, _dt: float, dx: float, Vp_sed, u_max: float, qs0: float,
         self._beta = beta
         self.stepmax = stepmax
         self.theta_sed = theta_sed
+        self.mod_erosion = mod_erosion
 
     def run(self, start_indices: np.ndarray, eta, stage, depth, cell_type,
             uw, ux, uy, Vp_dep_mud, Vp_dep_sand,
-            qw, qx, qy, qs) -> None:
+            qw, qx, qy, qs, mod_sed_weight) -> None:
         """The main function to route and deposit/erode sand parcels.
 
         Algorithm is to:
@@ -571,6 +583,7 @@ def run(self, start_indices: np.ndarray, eta, stage, depth, cell_type,
         self.pad_stage = shared_tools.custom_pad(stage)
         self.pad_depth = shared_tools.custom_pad(depth)
         self.pad_cell_type = shared_tools.custom_pad(cell_type)
+        self.pad_mod_sed_weight = shared_tools.custom_pad(mod_sed_weight)
         self.Vp_dep_mud = Vp_dep_mud
         self.Vp_dep_sand = Vp_dep_sand
         self.qw = qw
@@ -668,6 +681,7 @@ def _deposit_or_erode(self, px: int, py: int) -> None:
         qs_cap = (self.qs0 * self._f_bedload / self._u0**self._beta *
                   U_loc**self._beta)
         qs_loc = self.qs[px, py]
+        ero_mod_loc = self.mod_erosion[px, py]
 
         Vp_change = 0
         if qs_loc > qs_cap:
@@ -684,7 +698,8 @@ def _deposit_or_erode(self, px: int, py: int) -> None:
             #     critical erosion threshold for sand, *and* if the local
             #     transport capacity is not yet reached.
             Vp_change = self._compute_Vp_ero(self.Vp_sed, U_loc,
-                                             self.U_ero_sand, self._beta)
+                                             self.U_ero_sand, self._beta,
+                                             ero_mod_loc)
             Vp_change = self._limit_Vp_change(Vp_change, self.stage[px, py],
                                               self.eta[px, py], self._dx, 1)
             Vp_change = Vp_change * -1
@@ -716,7 +731,8 @@ class MudRouter(BaseRouter):
     """
     def __init__(self, _dt: float, dx: float, Vp_sed, u_max: float, U_dep_mud: float, U_ero_mud: float,
                  ivec_flat, jvec_flat, iwalk_flat, jwalk_flat, distances_flat,
-                 dry_depth: float, _lambda, beta: float, stepmax, theta_sed: float) -> None:
+                 dry_depth: float, _lambda, beta: float, stepmax,
+                 theta_sed: float, mod_erosion) -> None:
 
         self._dt = _dt
         self._dx = dx
@@ -735,10 +751,11 @@ def __init__(self, _dt: float, dx: float, Vp_sed, u_max: float, U_dep_mud: float
         self._beta = beta
         self.stepmax = stepmax
         self.theta_sed = theta_sed
+        self.mod_erosion = mod_erosion
 
     def run(self, start_indices, eta, stage, depth, cell_type,
             uw, ux, uy, Vp_dep_mud, Vp_dep_sand,
-            qw, qx, qy) -> None:
+            qw, qx, qy, mod_sed_weight) -> None:
         """The main function to route and deposit/erode mud parcels.
 
         """
@@ -753,6 +770,7 @@ def run(self, start_indices, eta, stage, depth, cell_type,
         self.pad_stage = shared_tools.custom_pad(stage)
         self.pad_depth = shared_tools.custom_pad(depth)
         self.pad_cell_type = shared_tools.custom_pad(cell_type)
+        self.pad_mod_sed_weight = shared_tools.custom_pad(mod_sed_weight)
         self.Vp_dep_mud = Vp_dep_mud
         self.Vp_dep_sand = Vp_dep_sand
         self.qw = qw
@@ -800,6 +818,7 @@ def _deposit_or_erode(self, px: int, py: int) -> None:
         .. important:: TODO: complete description specific for mud transport
         """
         U_loc = self.uw[px, py]
+        ero_mod_loc = self.mod_erosion[px, py]
 
         Vp_change = 0
         if U_loc < self.U_dep_mud:
@@ -811,7 +830,8 @@ def _deposit_or_erode(self, px: int, py: int) -> None:
 
         if U_loc > self.U_ero_mud:
             Vp_change = self._compute_Vp_ero(self.Vp_sed, U_loc,
-                                             self.U_ero_mud, self._beta)
+                                             self.U_ero_mud, self._beta,
+                                             ero_mod_loc)
             Vp_change = self._limit_Vp_change(Vp_change, self.stage[px, py],
                                               self.eta[px, py], self._dx, 1)
             Vp_change = Vp_change * -1