diff --git a/news/67.bugfix b/news/67.bugfix
new file mode 100644
index 0000000..82f64f4
--- /dev/null
+++ b/news/67.bugfix
@@ -0,0 +1,3 @@
+
+Fixed a bug where the topographic elevations were not being updated correctly after the
+components were time stepped.
diff --git a/news/67.feature b/news/67.feature
new file mode 100644
index 0000000..6b06227
--- /dev/null
+++ b/news/67.feature
@@ -0,0 +1,4 @@
+
+Added a :meth:`~sequence.SequenceModelGrid.get_profile` method to :class:`~sequence.SequenceModelGrid` that is a convenience
+function to return a field just along the middle row of nodes (i.e. along the
+profile).
diff --git a/news/67.feature.1 b/news/67.feature.1
new file mode 100644
index 0000000..01193e8
--- /dev/null
+++ b/news/67.feature.1
@@ -0,0 +1,6 @@
+
+Changed all components to only update the current amount of subsidence and
+deposited sediment but not to update secondary fields like the topographic
+elevation. This removes the dependence of these secondary fields on
+component ordering and ensures these secondary fields are updated once, and only
+once, with each time step.
diff --git a/news/67.misc b/news/67.misc
new file mode 100644
index 0000000..42332ce
--- /dev/null
+++ b/news/67.misc
@@ -0,0 +1,2 @@
+
+Added tests for the compaction, flexure, and subsidence components.
diff --git a/sequence/_grid.py b/sequence/_grid.py
index c944173..a422259 100644
--- a/sequence/_grid.py
+++ b/sequence/_grid.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 from landlab import RasterModelGrid
+from numpy.typing import NDArray
 
 if sys.version_info >= (3, 11):
     import tomllib
@@ -41,6 +42,21 @@ def __init__(self, n_cols: int, spacing: float = 100.0):
         self.at_node["sediment_deposit__thickness"] = self.zeros(at="node")
         self.at_grid["sea_level__elevation"] = 0.0
 
+    def get_profile(self, name: str) -> NDArray:
+        """Return the values of a field along the grid's profile.
+
+        Parameters
+        ----------
+        name : str
+            The name of an at-node field.
+
+        Returns
+        -------
+        values : ndarray
+            The values of the field located at the middle row of nodes.
+        """
+        return self.at_node[name].reshape(self.shape)[1]
+
     @classmethod
     def from_toml(cls, filepath: os.PathLike[str]) -> "SequenceModelGrid":
         """Load a :class:`~SequenceModelGrid` from a *toml*-formatted file.
diff --git a/sequence/cli.py b/sequence/cli.py
index 34dbb58..8f736af 100644
--- a/sequence/cli.py
+++ b/sequence/cli.py
@@ -231,8 +231,12 @@ def run(ctx: Any, with_citations: bool, dry_run: bool) -> None:
     run_dir = pathlib.Path.cwd()
 
     times, names = _find_config_files(".")
+    if len(times) == 0:
+        err("[ERROR] unable to find a configuration file.")
+        raise click.Abort()
+
     if not silent:
-        out(f"config_files = [{', '.join(repr(name) for name in names)}]")
+        out(f"[INFO] config_files: {', '.join(repr(name) for name in names)}")
     params = TimeVaryingConfig.from_files(names, times=times)
 
     model_params = params.as_dict()
@@ -250,20 +254,16 @@ def run(ctx: Any, with_citations: bool, dry_run: bool) -> None:
     )
 
     if verbose or not silent:
-        out("enabled = [")
         for name in model.components:
-            out(f"  {name!r},")
-        out("]")
-        out("disabled = [")
+            out(f"[INFO] ✅ enabled: {name}")
         for name in set(SequenceModel.ALL_PROCESSES) - set(model.components):
-            out(f"  {name!r},")
-        out("]")
+            out(f"[INFO] ❌ disabled {name}")
 
     if not silent and verbose:
         out(params.dump())
 
     if not silent and len(processes) == 0:
-        out("⚠️  ALL PROCESSES HAVE BEEN DISABLED! ⚠️")
+        out("[WARN] ⚠️  ALL PROCESSES HAVE BEEN DISABLED! ⚠️")
 
     if not silent and with_citations:
         from landlab.core.model_component import registry
diff --git a/sequence/fluvial.py b/sequence/fluvial.py
index d258359..ebaf6ea 100644
--- a/sequence/fluvial.py
+++ b/sequence/fluvial.py
@@ -34,6 +34,14 @@ class Fluvial(Component):
             "mapping": "grid",
             "doc": "Position of sea level",
         },
+        "sediment_deposit__thickness": {
+            "dtype": "float",
+            "intent": "inout",
+            "optional": False,
+            "units": "m",
+            "mapping": "node",
+            "doc": "Thickness of deposition or erosion",
+        },
         "delta_sediment_sand__volume_fraction": {
             "dtype": "float",
             "intent": "out",
@@ -135,17 +143,17 @@ def run_one_step(self, dt: float) -> None:
         channel_width = sand_vol / qs
 
         x = self.grid.x_of_node.reshape(self.grid.shape)[1]
-        z = self.grid.at_node["topographic__elevation"].reshape(self.grid.shape)[1]
+        elevation = self.grid.get_profile("topographic__elevation")
 
         shore = find_shoreline(
-            x, z, sea_level=self.grid.at_grid["sea_level__elevation"]
+            x, elevation, sea_level=self.grid.at_grid["sea_level__elevation"]
         )
 
         land = x < shore
         water = x >= shore
         # land = self.grid.x_of_node[self.grid.node_at_cell] < shore
         # slope = np.gradient(z[1, land]) / self.grid.dx
-        slope = np.gradient(z) / self.grid.dx
+        slope = np.gradient(elevation) / self.grid.dx
         # slp = np.zeros(1)
         # slp[0] = self.plain_slope
         # slope = np.concatenate((slp,slop))
@@ -155,7 +163,7 @@ def run_one_step(self, dt: float) -> None:
         # )
 
         # channel_depth = np.zeros(self.grid.number_of_cells)
-        channel_depth = np.zeros_like(z)  # use upsteam slope
+        channel_depth = np.zeros_like(elevation)  # use upsteam slope
         # channel_depth[0] = (
         #    (self.sand_density - 1000.) / 1000. * self.sand_grain / self.plain_slope
         channel_depth[land] = (
@@ -170,28 +178,13 @@ def run_one_step(self, dt: float) -> None:
         # width_cb = channel_width/epsilon
 
         # Original: r_cb = (model.new_height[i]-model.height[i]+model.d_sl);
-        r_cb = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-            * epsilon
-        )
-        dz = (
-            self.grid.at_node["bedrock_surface__increment_of_elevation"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
+        r_cb = self.grid.get_profile("sediment_deposit__thickness") * epsilon
+        dz = self.grid.get_profile("bedrock_surface__increment_of_elevation").copy()
         # original: r_b = model.thickness[i];
-        r_b = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
+        r_b = self.grid.get_profile("sediment_deposit__thickness").copy()
         # r_fp = np.zeros(self.grid.shape[1])
-        r_fp = np.zeros_like(z)
-        percent_sand = self.grid.at_node[
-            "delta_sediment_sand__volume_fraction"
-        ].reshape(self.grid.shape)[1]
+        r_fp = np.zeros_like(elevation)
+        percent_sand = self.grid.get_profile("delta_sediment_sand__volume_fraction")
         percent_sand.fill(0.0)
 
         for i in np.where(land)[0]:
@@ -295,25 +288,13 @@ def run_one_step(self, dt: float) -> None:
                 if add_mud[i] < 0.0:
                     add_mud[i] = 0.0
 
-        thickness = (
-            self.grid.at_node["sediment_deposit__thickness"]
-            .reshape(self.grid.shape)[1]
-            .copy()
-        )
-        thickness[water] += add_mud[water]
+        sediment_thickness = self.grid.get_profile("sediment_deposit__thickness")
+        sediment_thickness[water] += add_mud[water]
 
         np.divide(
-            thickness[water] - add_mud[water],
-            thickness[water],
-            where=thickness[water] > 0.0,
+            sediment_thickness[water] - add_mud[water],
+            sediment_thickness[water],
+            where=sediment_thickness[water] > 0.0,
             out=percent_sand[water],
         )
-
         np.clip(percent_sand[water], 0.0, 1.0, out=percent_sand[water])
-
-        plus_mud = np.zeros(self.grid.shape)
-        plus_mud[1] = add_mud
-        sdt = self.grid.at_node["sediment_deposit__thickness"].reshape(plus_mud.shape)
-        sdt[1][:] += plus_mud[1]
-        te = self.grid.at_node["topographic__elevation"].reshape(plus_mud.shape)
-        te[1][:] -= plus_mud[1]
diff --git a/sequence/sediment_flexure.py b/sequence/sediment_flexure.py
index a1ebe65..b207700 100644
--- a/sequence/sediment_flexure.py
+++ b/sequence/sediment_flexure.py
@@ -1,6 +1,9 @@
 """Subside a `SequenceModelGrid` using flexure."""
+from typing import Union
+
 import numpy as np
 from landlab.components.flexure import Flexure1D
+from numpy.typing import NDArray
 
 from ._grid import SequenceModelGrid
 
@@ -33,35 +36,27 @@ class SedimentFlexure(Flexure1D):
         },
         "topographic__elevation": {
             "dtype": "float",
-            "intent": "inout",
+            "intent": "in",
             "optional": True,
             "units": "m",
             "mapping": "node",
             "doc": "land and ocean bottom elevation, positive up",
         },
-        "bedrock_surface__elevation": {
+        "delta_sediment_sand__volume_fraction": {
             "dtype": "float",
-            "intent": "inout",
+            "intent": "in",
             "optional": True,
-            "units": "m",
+            "units": "-",
             "mapping": "node",
-            "doc": "New bedrock elevation following subsidence",
+            "doc": "delta sand fraction",
         },
         "bedrock_surface__increment_of_elevation": {
             "dtype": "float",
-            "intent": "out",
+            "intent": "inout",
             "optional": True,
             "units": "m",
             "mapping": "node",
-            "doc": "Amount of subsidence due to sediment loading",
-        },
-        "delta_sediment_sand__volume_fraction": {
-            "dtype": "float",
-            "intent": "in",
-            "optional": True,
-            "units": "-",
-            "mapping": "node",
-            "doc": "delta sand fraction",
+            "doc": "Total amount of subsidence",
         },
     }
 
@@ -70,7 +65,7 @@ def __init__(
         grid: SequenceModelGrid,
         sand_density: float = 2650,
         mud_density: float = 2720.0,
-        isostasytime: float = 7000.0,
+        isostasytime: Union[float, None] = 7000.0,
         water_density: float = 1030.0,
         **kwds: dict,
         # **sediments,
@@ -101,13 +96,16 @@ def __init__(
             self.mud_density, self.water_density, 0.65
         )
 
-        self._isostasytime = SedimentFlexure.validate_isostasy_time(isostasytime)
-
-        self._dt = 100.0  # default timestep = 100 y
-
-        # grid.add_zeros("lithosphere__increment_of_overlying_pressure", at="node")
+        if np.isclose(isostasytime, 0.0):
+            isostasytime = None
+        self._isostasy_time = (
+            SedimentFlexure.validate_isostasy_time(isostasytime)
+            if isostasytime is not None
+            else None
+        )
+        self._dt = 1.0
 
-        super().__init__(grid, **kwds)
+        super().__init__(grid, rows=1, **kwds)
 
         if "lithosphere__increment_of_overlying_pressure" not in grid.at_node:
             grid.add_zeros("lithosphere__increment_of_overlying_pressure", at="node")
@@ -120,7 +118,9 @@ def __init__(
         if "bedrock_surface__increment_of_elevation" not in grid.at_node:
             grid.add_empty("bedrock_surface__increment_of_elevation", at="node")
 
-        self.subs_pool = self.grid.zeros(at="node")
+        self.subs_pool = np.zeros_like(
+            self.grid.get_profile("lithosphere_surface__increment_of_elevation")
+        )
 
     @staticmethod
     def _calc_bulk_density(
@@ -174,6 +174,11 @@ def validate_isostasy_time(time: float) -> float:
             raise ValueError(f"negative isostasy time ({time})")
         return time
 
+    @property
+    def isostasy_time(self) -> Union[float, None]:
+        """Return the isostasy time."""
+        return self._isostasy_time
+
     @property
     def sand_density(self) -> float:
         """Return the density of sand."""
@@ -194,7 +199,7 @@ def sand_bulk_density(self) -> float:
 
     @property
     def mud_density(self) -> float:
-        """Return teh density of mud."""
+        """Return the density of mud."""
         return self._mud_density
 
     @mud_density.setter
@@ -225,50 +230,144 @@ def water_density(self, density: float) -> None:
             self.mud_density, self.water_density, 0.65
         )
 
-    def update(self) -> None:
-        """Update the component by a single time step."""
-        if self._isostasytime > 0.0:
-            isostasyfrac = 1 - np.exp(-1.0 * self._dt / self._isostasytime)
-        else:
-            isostasyfrac = 1.0
+    @staticmethod
+    def _calc_loading(
+        deposit_thickness: NDArray[np.floating],
+        z: NDArray[np.floating],
+        sediment_density: Union[float, NDArray[np.floating]],
+        water_density: float,
+    ) -> NDArray[np.floating]:
+        """Calculate the loading due to a, possibly submerged, sediment deposit.
+
+        Parameters
+        ----------
+        deposit_thickness : array-like
+            The amount of sediment deposited along the profile.
+        z : array-like
+            The elevation of the profile before the new sediment has
+            been added.
+        sediment_density : float
+            The bulk density of the added sediment.
+        water_density : float
+            The density of water.
 
-        self.grid.at_node["lithosphere_surface__increment_of_elevation"][:] = 0.0
+        Returns
+        -------
+        ndarray
+            The loading along the profile.
+        """
+        z_new = z + deposit_thickness
 
-        # calculate density based on sand_frac
-        percent_sand = self.grid.at_node[
-            "delta_sediment_sand__volume_fraction"
-        ].reshape(self.grid.shape)[1]
-        rho_sediment = (
-            percent_sand * self._rho_sand + (1 - percent_sand) * self._rho_mud
+        dry = (z >= 0.0) & (z_new >= 0.0)
+        wet = (z <= 0.0) & (z_new <= 0.0)
+        mixed = ~dry & ~wet
+
+        density = np.full_like(z, sediment_density)
+        density[wet] -= water_density
+
+        dry_fraction = np.abs(
+            np.maximum(z_new[mixed], z[mixed]) / deposit_thickness[mixed]
         )
+        wet_fraction = 1.0 - dry_fraction
+
+        density[mixed] -= wet_fraction * water_density
+
+        return density * deposit_thickness
+
+    @staticmethod
+    def _calc_density(
+        sand_fraction: Union[float, NDArray], sand_density: float, mud_density: float
+    ) -> Union[float, NDArray]:
+        """Calculate the density of a sand/mud mixture.
+
+        Parameters
+        ----------
+        sand_fraction : float or ndarray
+            Fraction of sediment that is sand. The rest is mud.
+        sand_density : float
+            The density of the sand.
+        mud_density : float
+            The density of the mud.
+
+        Returns
+        -------
+        float or ndarray
+            The density of the mixture.
+
+        Examples
+        --------
+        >>> from sequence.sediment_flexure import SedimentFlexure
+        >>> SedimentFlexure._calc_density(0.5, 1600, 1200)
+        1400.0
+        >>> SedimentFlexure._calc_density([1.0, 0.75, 0.5, 0.25, 0.0], 1600.0, 1200.0)
+        array([ 1600.,  1500.,  1400.,  1300.,  1200.])
+        """
+        sand_fraction = np.asarray(sand_fraction)
+        return sand_fraction * sand_density + (1.0 - sand_fraction) * mud_density
+
+    @staticmethod
+    def _calc_isostasy_fraction(isostasy_time: Union[float, None], dt: float) -> float:
+        """Calculate the fraction of isostatic subsidence.
+
+        Parameters
+        ----------
+        isostasy_time : float or None
+            The time parameter that represente the e-folding time
+            to isostasy. If ``None``, assume isostasy.
+        dt : float
+            Time step.
+
+        Returns
+        -------
+        float
+            Fraction of isostatic subsidence that occurs over the time step.
+        """
+        if isostasy_time is None:
+            return 1.0
+        else:
+            return 1.0 - np.exp(-dt / isostasy_time)
+
+    def update(self) -> None:
+        """Update the component by a single time step."""
+        self.grid.get_profile("lithosphere_surface__increment_of_elevation").fill(0.0)
 
-        # load underwater displaces water
-        z = self.grid.at_node["topographic__elevation"].reshape(self.grid.shape)[1]
         sea_level = self.grid.at_grid["sea_level__elevation"]
-        under_water = z < sea_level
-        rho_sediment[under_water] = rho_sediment[under_water] - self.water_density
+        elevation = self.grid.get_profile("topographic__elevation")
+        deposit_thickness = self.grid.get_profile("sediment_deposit__thickness")
+        sand_fraction = self.grid.get_profile("delta_sediment_sand__volume_fraction")
 
-        dz = self.grid.at_node["sediment_deposit__thickness"].reshape(self.grid.shape)[
-            1
-        ]
-        pressure = self.grid.at_node[
-            "lithosphere__increment_of_overlying_pressure"
-        ].reshape(self.grid.shape)[1]
-        pressure[:] = dz * rho_sediment * self.gravity * self.grid.dx
+        sediment_density = self._calc_density(
+            sand_fraction, self.sand_bulk_density, self.mud_bulk_density
+        )
 
-        Flexure1D.update(self)
+        pressure = self.grid.get_profile("lithosphere__increment_of_overlying_pressure")
+        pressure[:] = (
+            self._calc_loading(
+                deposit_thickness,
+                elevation - sea_level,
+                sediment_density,
+                self.water_density,
+            )
+            * self.gravity
+            * self.grid.dx
+        )
 
-        dz = self.grid.at_node["lithosphere_surface__increment_of_elevation"]
-        self.subs_pool[:] += dz
-        dz = self.subs_pool[:] * isostasyfrac
-        self.subs_pool[:] = self.subs_pool[:] - dz
+        Flexure1D.update(self)
 
-        self.grid.at_node["bedrock_surface__increment_of_elevation"][:] = dz
+        isostatic_deflection = self.grid.get_profile(
+            "lithosphere_surface__increment_of_elevation"
+        )
+        isostasy_fraction = self._calc_isostasy_fraction(self.isostasy_time, self._dt)
+        self.subs_pool[:] += isostatic_deflection
+        deflection = self.subs_pool[:] * isostasy_fraction
+        self.subs_pool[:] = self.subs_pool[:] - deflection
 
-        self.grid.at_node["bedrock_surface__elevation"] -= dz
-        self.grid.at_node["topographic__elevation"] -= dz
+        total_deflection = self.grid.get_profile(
+            "bedrock_surface__increment_of_elevation"
+        )
+        total_deflection += deflection
 
-    def run_one_step(self, dt: float = 100.0) -> None:
+    def run_one_step(self, dt: float = 1.0) -> None:
         """Update the component by a time step.
 
         Parameters
diff --git a/sequence/sequence.py b/sequence/sequence.py
index a5e5f88..38825b6 100644
--- a/sequence/sequence.py
+++ b/sequence/sequence.py
@@ -90,7 +90,14 @@ def __init__(
             - self.grid.at_node["bedrock_surface__elevation"]
         )
 
-        self.add_layer(initial_sediment_thickness[self.grid.node_at_cell])
+        if (initial_sediment_thickness[self.grid.node_at_cell] > 0.0).any():
+            self.add_layer(initial_sediment_thickness[self.grid.node_at_cell])
+
+        self.grid.at_node["topographic__elevation"][self.grid.node_at_cell] = (
+            self.grid.at_node["bedrock_surface__elevation"][self.grid.node_at_cell]
+            + self.grid.event_layers.thickness
+        )
+        self.grid.at_node["sediment_deposit__thickness"].fill(0.0)
 
     @property
     def time(self) -> float:
@@ -116,6 +123,11 @@ def update(self, dt: Optional[float] = None) -> None:
         self.add_layer(
             self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell]
         )
+        self.grid.at_node["topographic__elevation"][self.grid.node_at_cell] = (
+            self.grid.at_node["bedrock_surface__elevation"][self.grid.node_at_cell]
+            + self.grid.event_layers.thickness
+        )
+        self.grid.at_node["sediment_deposit__thickness"].fill(0.0)
 
     def run(
         self,
@@ -214,6 +226,10 @@ def add_layer(self, dz_at_cell: NDArray[np.floating]) -> None:
         except KeyError:
             x_of_shelf_edge = np.nan
 
+        self.grid.at_node["topographic__elevation"][
+            self.grid.node_at_cell
+        ] += dz_at_cell
+
         self.grid.event_layers.add(dz_at_cell, **self.layer_properties())
         self.grid.at_layer_grid.add(
             1.0,
diff --git a/sequence/sequence_model.py b/sequence/sequence_model.py
index 4a604bf..52e0cf2 100644
--- a/sequence/sequence_model.py
+++ b/sequence/sequence_model.py
@@ -2,12 +2,12 @@
 import warnings
 from collections import OrderedDict, defaultdict
 from collections.abc import Hashable, Iterable
-from typing import Dict, Optional
+from typing import Any, Dict, Optional
 
 import numpy as np
 from compaction.landlab import Compact
-from landlab import FieldError
 from landlab.bmi.bmi_bridge import TimeStepper
+from numpy.typing import ArrayLike
 
 from ._grid import SequenceModelGrid
 from .bathymetry import BathymetryReader
@@ -256,46 +256,95 @@ def advance_components(self, dt: float) -> None:
         dt : float
             The time step to advance the components.
         """
+        self.grid.at_node["sediment_deposit__thickness"].fill(0.0)
+
         for component in self._components.values():
             component.run_one_step(dt)
 
+        self._update_fields()
+
+        self.grid.event_layers.add(
+            self.grid.at_node["sediment_deposit__thickness"][self.grid.node_at_cell],
+            **self.layer_properties(),
+        )
+
+        if (
+            self.grid.event_layers.number_of_layers - self._n_archived_layers
+        ) % 20 == 0:
+            self.grid.event_layers.reduce(
+                self._n_archived_layers,
+                self._n_archived_layers + 10,
+                **self.layer_reducers(),
+            )
+            self._n_archived_layers += 1
+
+    def layer_properties(self) -> dict[str, ArrayLike]:
+        """Return the properties being tracked at each layer.
+
+        Returns
+        -------
+        properties : dict
+            A dictionary of the tracked properties where the keys
+            are the names of properties and the values are the
+            property values at each column.
+        """
         dz = self.grid.at_node["sediment_deposit__thickness"]
-        # percent_sand = self.grid.at_node["delta_sediment_sand__volume_fraction"]
         water_depth = (
             self.grid.at_grid["sea_level__elevation"]
             - self.grid.at_node["topographic__elevation"]
         )
 
-        layer_properties = dict(
-            age=self.clock.time,
-            water_depth=water_depth[self.grid.node_at_cell],
-            t0=dz[self.grid.node_at_cell].clip(0.0),
-            # percent_sand=percent_sand[self.grid.node_at_cell],
-        )
+        properties = {
+            "age": self.clock.time,
+            "water_depth": water_depth[self.grid.node_at_cell],
+            "t0": dz[self.grid.node_at_cell].clip(0.0),
+            "porosity": 0.5,
+        }
+
         if "compaction" in self._components:
-            layer_properties["porosity"] = self._components["compaction"].porosity_max
+            properties["porosity"] = self._components["compaction"].porosity_max
+
         try:
             percent_sand = self.grid.at_node["delta_sediment_sand__volume_fraction"]
-        except FieldError:
+        except KeyError:
             pass
         else:
-            layer_properties["percent_sand"] = percent_sand[self.grid.node_at_cell]
+            properties["percent_sand"] = percent_sand[self.grid.node_at_cell]
 
-        self.grid.event_layers.add(dz[self.grid.node_at_cell], **layer_properties)
+        return properties
 
-        if (
-            self.grid.event_layers.number_of_layers - self._n_archived_layers
-        ) % 20 == 0:
-            self.grid.event_layers.reduce(
-                self._n_archived_layers,
-                self._n_archived_layers + 10,
-                age=np.max,
-                percent_sand=np.mean,
-                porosity=np.mean,
-                t0=np.sum,
-                water_depth=np.mean,
-            )
-            self._n_archived_layers += 1
+    def layer_reducers(self) -> dict[str, Any]:
+        """Return layer-reducers for each property.
+
+        Returns
+        -------
+        reducers : dict
+            A dictionary of reducers where keys are property names and values
+            are functions that act as layer reducers for those quantities.
+        """
+        reducers = {
+            "age": np.max,
+            "water_depth": np.mean,
+            "t0": np.sum,
+            "porosity": np.mean,
+        }
+        if "percent_sand" in self.grid.event_layers.tracking:
+            reducers["percent_sand"] = np.mean
+
+        return reducers
+
+    def _update_fields(self) -> None:
+        """Update fields that depend on other fields."""
+        if "bedrock_surface__increment_of_elevation" in self.grid.at_node:
+            self.grid.at_node["bedrock_surface__elevation"] += self.grid.at_node[
+                "bedrock_surface__increment_of_elevation"
+            ]
+            self.grid.at_node["bedrock_surface__increment_of_elevation"][:] = 0.0
+
+        self.grid.at_node["topographic__elevation"][self.grid.node_at_cell] = (
+            self.grid.at_node["bedrock_surface__elevation"][self.grid.node_at_cell]
+            + self.grid.event_layers.thickness
+        )
 
 
 def _match_values(d1: dict, d2: dict, keys: Iterable[Hashable]) -> None:
diff --git a/sequence/submarine.py b/sequence/submarine.py
index 68b73c1..4284b38 100644
--- a/sequence/submarine.py
+++ b/sequence/submarine.py
@@ -35,7 +35,7 @@ class SubmarineDiffuser(LinearDiffuser):
         },
         "sediment_deposit__thickness": {
             "dtype": "float",
-            "intent": "out",
+            "intent": "inout",
             "optional": False,
             "units": "m",
             "mapping": "node",
@@ -271,11 +271,9 @@ def run_one_step(self, dt: float) -> None:
         dt : float
             Time step to advance component by.
         """
-        z_before = self.grid.at_node["topographic__elevation"].copy()
-
         shore = find_shoreline(
             self.grid.x_of_node[self.grid.node_at_cell],
-            z_before[self.grid.node_at_cell],
+            self.grid.at_node["topographic__elevation"][self.grid.node_at_cell],
             sea_level=self.grid.at_grid["sea_level__elevation"],
         )
 
@@ -289,10 +287,13 @@ def run_one_step(self, dt: float) -> None:
         z[1, 0] = z[1, 1] + self._plain_slope * (x[1, 1] - x[1, 0])
         # self._load/self._load0)
 
+        z_before = self.grid.at_node["topographic__elevation"].copy()
+
         super().run_one_step(dt)
 
-        self.grid.at_node["sediment_deposit__thickness"][:] = (
-            self.grid.at_node["topographic__elevation"] - z_before
-        )
+        dz = self.grid.at_node["topographic__elevation"] - z_before
+        self.grid.at_node["topographic__elevation"][:] = z_before
+
+        self.grid.at_node["sediment_deposit__thickness"][:] += dz
 
         self._time += dt
diff --git a/sequence/subsidence.py b/sequence/subsidence.py
index b0e131a..e4c949f 100644
--- a/sequence/subsidence.py
+++ b/sequence/subsidence.py
@@ -19,21 +19,13 @@ class SubsidenceTimeSeries(Component):
 
     _info = {
         "bedrock_surface__increment_of_elevation": {
-            "dtype": "float",
-            "intent": "out",
-            "optional": False,
-            "units": "m",
-            "mapping": "node",
-            "doc": "Increment of elevation",
-        },
-        "bedrock_surface__elevation": {
             "dtype": "float",
             "intent": "inout",
             "optional": False,
             "units": "m",
             "mapping": "node",
-            "doc": "Surface elevation",
-        },
+            "doc": "Change in elevation due to subsidence",
+        }
     }
 
     def __init__(
@@ -52,8 +44,11 @@ def __init__(
             'nearest', 'zero', 'slinear', 'quadratic', 'cubic').
             Default is 'linear'.
         """
-        if "bedrock_surface__elevation" not in grid.at_node:
-            grid.add_empty("bedrock_surface__elevation", at="node")
+        if "bedrock_surface__increment_of_elevation" not in grid.at_node:
+            grid.add_zeros("bedrock_surface__increment_of_elevation", at="node")
+
+        if "lithosphere_surface__increment_of_elevation" not in grid.at_node:
+            grid.add_zeros("lithosphere_surface__increment_of_elevation", at="node")
 
         super().__init__(grid)
 
@@ -61,20 +56,19 @@ def __init__(
         self._kind = kind
 
         data = np.loadtxt(filepath, delimiter=",", comments="#")
-        subsidence = SubsidenceTimeSeries._subsidence_interpolator(
+        self._subsidence = SubsidenceTimeSeries._subsidence_interpolator(
             data, kind=self._kind
         )
-        if "bedrock_surface__increment_of_elevation" not in grid.at_node:
-            self.grid.add_empty("bedrock_surface__increment_of_elevation", at="node")
 
-        inc = self.grid.at_node["bedrock_surface__increment_of_elevation"].reshape(
-            self.grid.shape
-        )
-        inc[:] = subsidence(self.grid.x_of_node[self.grid.nodes_at_bottom_edge])
+        self._dz_dt = self._calc_subsidence_rate()
 
-        self._dz = inc.copy()
         self._time = 0.0
 
+    @property
+    def subsidence_rate(self) -> NDArray:
+        """Return the current subsidence rate."""
+        return self._dz_dt
+
     @staticmethod
     def _subsidence_interpolator(
         data: NDArray, kind: str = "linear"
@@ -88,6 +82,9 @@ def _subsidence_interpolator(
             bounds_error=True,
         )
 
+    def _calc_subsidence_rate(self) -> NDArray:
+        return self._subsidence(self.grid.x_of_node[self.grid.nodes_at_bottom_edge])
+
     @property
     def time(self) -> float:
         """Return the current component time."""
@@ -101,14 +98,11 @@ def filepath(self) -> str:
     @filepath.setter
     def filepath(self, new_path: os.PathLike) -> None:
         self._filepath = new_path
-        subsidence = SubsidenceTimeSeries._subsidence_interpolator(
+        self._subsidence = SubsidenceTimeSeries._subsidence_interpolator(
             np.loadtxt(self._filepath, delimiter=",", comments="#"), kind=self._kind
         )
-        inc = self.grid.at_node["bedrock_surface__increment_of_elevation"].reshape(
-            self.grid.shape
-        )
-        inc[:] = subsidence(self.grid.x_of_node[self.grid.nodes_at_bottom_edge])
-        self._dz = inc.copy()
+
+        self._dz_dt = self._calc_subsidence_rate()
 
     def run_one_step(self, dt: float) -> None:
         """Update the component by a time step.
@@ -118,16 +112,8 @@ def run_one_step(self, dt: float) -> None:
         dt : float
             The time step to update the component by.
         """
-        dz = self.grid.at_node["bedrock_surface__increment_of_elevation"]
-        z = self.grid.at_node["bedrock_surface__elevation"]
-        z_top = self.grid.at_node["topographic__elevation"]
-
-        dz = dz.reshape(self.grid.shape)
-        z = z.reshape(self.grid.shape)
-        z_top = z_top.reshape(self.grid.shape)
-
-        dz[:] = self._dz * dt
-        z[:] += dz
-        z_top[:] += dz
+        self.grid.get_profile("bedrock_surface__increment_of_elevation")[:] += (
+            self.subsidence_rate * dt
+        )
 
         self._time += dt
diff --git a/tests/test_compaction.py b/tests/test_compaction.py
new file mode 100644
index 0000000..3ef31d5
--- /dev/null
+++ b/tests/test_compaction.py
@@ -0,0 +1,65 @@
+from numpy.testing import assert_array_almost_equal
+from pytest import approx
+
+from sequence import Sequence, SequenceModelGrid
+from sequence.processes import Compact
+
+
+def test_layers_compact():
+    grid = SequenceModelGrid(5)
+    grid.add_zeros("bedrock_surface__elevation", at="node")
+    grid.add_zeros("topographic__elevation", at="node")
+
+    compaction = Compact(
+        grid,
+        c=5e-08,
+        porosity_max=0.5,
+        porosity_min=0.01,
+        rho_grain=2650.0,
+        rho_void=1000.0,
+    )
+    sequence = Sequence(grid, time_step=100.0, components=[compaction])
+    assert grid.event_layers.number_of_layers == 0
+    for _ in range(n_layers := 15):
+        sequence.add_layer(100.0)
+    assert grid.event_layers.number_of_layers == 15
+    assert grid.at_node["topographic__elevation"][grid.node_at_cell] == approx(
+        n_layers * 100.0
+    )
+
+    sequence.run(progress_bar=False)
+
+    assert (grid.event_layers["porosity"][: n_layers - 1, :] < 0.5).all()
+    assert grid.event_layers["porosity"][n_layers - 1, :] == approx(0.5)
+
+    assert (grid.event_layers.dz[: n_layers - 1, :] < 100.0).all()
+    assert grid.event_layers.dz[n_layers - 1, :] == approx(100.0)
+
+
+def test_layers_compact_only_once():
+    grid = SequenceModelGrid(5)
+    grid.add_zeros("bedrock_surface__elevation", at="node")
+    grid.add_zeros("topographic__elevation", at="node")
+
+    compaction = Compact(
+        grid,
+        c=5e-08,
+        porosity_max=0.5,
+        porosity_min=0.01,
+        rho_grain=2650.0,
+        rho_void=1000.0,
+    )
+    sequence = Sequence(grid, time_step=100.0, components=[compaction])
+    assert grid.event_layers.number_of_layers == 0
+    for _ in range(15):
+        sequence.add_layer(100.0)
+
+    sequence.run(progress_bar=False)
+    assert sequence.time == approx(100.0)
+
+    thickness_of_layers = grid.event_layers.thickness.copy()
+
+    sequence.run(until=4000.0, progress_bar=False)
+    assert sequence.time == approx(4000.0)
+
+    assert_array_almost_equal(grid.event_layers.thickness, thickness_of_layers)
diff --git a/tests/test_flexure.py b/tests/test_flexure.py
index 47be240..2b999e4 100644
--- a/tests/test_flexure.py
+++ b/tests/test_flexure.py
@@ -1,15 +1,15 @@
+import numpy as np
 import pytest
-from landlab import RasterModelGrid
+from numpy.testing import assert_array_almost_equal
 
+from sequence import SequenceModelGrid
 from sequence.sediment_flexure import SedimentFlexure
 
 
 @pytest.fixture()
 def grid():
-    grid = RasterModelGrid((3, 4))
-    grid.add_empty("sediment_deposit__thickness", at="node")
+    grid = SequenceModelGrid(4, spacing=1.0)
     grid.add_empty("topographic__elevation", at="node")
-    grid.add_empty("bedrock_surface__elevation", at="node")
     return grid
 
 
@@ -61,18 +61,91 @@ def test_update_bulk_density(grid, prop):
 
 
 def test_flexure():
-    grid = RasterModelGrid((3, 4))
-    grid.add_empty("sediment_deposit__thickness", at="node")
-    grid.add_zeros("topographic__elevation", at="node")
-    grid.add_zeros("bedrock_surface__elevation", at="node")
+    grid = SequenceModelGrid(100)
+    initial_elevation = grid.add_zeros("topographic__elevation", at="node").copy()
 
     grid.at_node["sediment_deposit__thickness"][:] = 1.0
     flexure = SedimentFlexure(grid)
     flexure.run_one_step()
 
-    assert (
-        grid.at_node["topographic__elevation"].reshape(grid.shape)[1, 1:-1] < 0.0
-    ).all()
-    assert (
-        grid.at_node["bedrock_surface__elevation"].reshape(grid.shape)[1, 1:-1] < 0.0
-    ).all()
+    assert (grid.get_profile("lithosphere_surface__increment_of_elevation") > 0.0).all()
+    assert (grid.at_node["topographic__elevation"] == initial_elevation).all()
+
+
+def test_all_dry():
+    z = np.asarray([4.0, 3.0, 2.0, 1.0, 0.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.full_like(z, 2.0)
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, dz * sediment_density)
+
+
+def test_all_wet():
+    z = np.asarray([-2.0, -3.0, -4.0, -5.0, -6.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.full_like(z, 2.0)
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, dz * water_density)
+
+
+def test_loading_wet_or_dry():
+    z = np.asarray([2.0, 1.0, 0.0, -1.0, -2.0])
+    sediment_density = 2000.0
+    water_density = 1000.0
+    dz = np.full_like(z, 0.5)
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(
+        loading, [0.5 * 2000.0, 0.5 * 2000, 0.5 * 2000, 0.5 * 1000, 0.5 * 1000]
+    )
+
+
+def test_loading_mixed():
+    z = np.asarray([1.0, 0.0, -1.0, -2.0, -3.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.full_like(z, 2.0)
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, [2.0, 2.0, 1.0 + 0.5, 1.0, 1.0])
+
+
+def test_dry_erosion():
+    z = np.asarray([4.0, 3.0, 2.0, 1.0, 0.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.asarray([-1.0, -1.0, 1.0, 1.0, 1.0])
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, dz * sediment_density)
+
+
+def test_wet_erosion():
+    z = np.asarray([0.0, -1.0, -2.0, -3.0, -4.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.asarray([-1.0, -1.0, 1.0, 1.0, 1.0])
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, dz * (sediment_density - water_density))
+
+
+def test_mixed_erosion():
+    z = np.asarray([3.0, 2.0, 1.0, 0.0, -1.0])
+    sediment_density = 1.0
+    water_density = 0.5
+    dz = np.asarray([-2.0, -2.0, -2.0, 2.0, 2.0])
+
+    loading = SedimentFlexure._calc_loading(dz, z, sediment_density, water_density)
+
+    assert_array_almost_equal(loading, [-2.0, -2.0, -1.0 - 0.5, 2.0, 1.0 + 0.5])
diff --git a/tests/test_subsidence.py b/tests/test_subsidence.py
index 4c44c2c..180be0b 100644
--- a/tests/test_subsidence.py
+++ b/tests/test_subsidence.py
@@ -1,14 +1,12 @@
 import pytest
-from landlab import RasterModelGrid
 from numpy.testing import assert_array_equal
 
+from sequence import SequenceModelGrid
 from sequence.subsidence import SubsidenceTimeSeries
 
 
 def test_bad_filepath(tmpdir):
-    grid = RasterModelGrid((3, 5))
-    grid.add_full("bedrock_surface__elevation", 0.0, at="node")
-    grid.add_full("topographic__elevation", 0.0, at="node")
+    grid = SequenceModelGrid(5)
 
     with pytest.raises(TypeError):
         SubsidenceTimeSeries(grid)
@@ -18,85 +16,108 @@ def test_bad_filepath(tmpdir):
             SubsidenceTimeSeries(grid, filepath="missing-file.csv")
 
 
+def test_time_step(tmpdir):
+    grid = SequenceModelGrid(5, spacing=1.0)
+
+    with tmpdir.as_cwd():
+        with open("subsidence.csv", "w") as fp:
+            fp.write("0.0,1.0\n5.0,1.0")
+        subsidence = SubsidenceTimeSeries(grid, filepath="subsidence.csv")
+        assert_array_equal(subsidence.subsidence_rate, 1.0)
+        assert_array_equal(grid.at_node["bedrock_surface__increment_of_elevation"], 0.0)
+
+        subsidence.run_one_step(dt=1.0)
+        assert_array_equal(subsidence.subsidence_rate, 1.0)
+        assert_array_equal(
+            grid.get_profile("bedrock_surface__increment_of_elevation"), 1.0
+        )
+
+        grid.at_node["bedrock_surface__increment_of_elevation"][:] = 0.0
+        subsidence.run_one_step(dt=2.0)
+        assert_array_equal(
+            grid.get_profile("bedrock_surface__increment_of_elevation"), 2.0
+        )
+
+
 def test_constant_subsidence(tmpdir):
-    grid = RasterModelGrid((3, 5))
-    grid.add_full("bedrock_surface__elevation", 0.0, at="node")
-    grid.add_full("topographic__elevation", 0.0, at="node")
+    grid = SequenceModelGrid(5, spacing=1.0)
 
     with tmpdir.as_cwd():
         with open("subsidence.csv", "w") as fp:
             fp.write("0.0,1.0\n5.0,1.0")
         subsidence = SubsidenceTimeSeries(grid, filepath="subsidence.csv")
-        assert_array_equal(grid.at_node["bedrock_surface__increment_of_elevation"], 1.0)
+        assert_array_equal(grid.at_node["bedrock_surface__increment_of_elevation"], 0.0)
 
         subsidence.run_one_step(dt=1.0)
-        assert_array_equal(grid.at_node["bedrock_surface__elevation"], 1.0)
+        assert_array_equal(
+            grid.get_profile("bedrock_surface__increment_of_elevation"), 1.0
+        )
 
+        grid.at_node["bedrock_surface__increment_of_elevation"][:] = 0.0
         subsidence.run_one_step(dt=1.0)
-        assert_array_equal(grid.at_node["bedrock_surface__elevation"], 2.0)
+        assert_array_equal(
+            grid.get_profile("bedrock_surface__increment_of_elevation"), 1.0
+        )
 
 
 def test_linear_subsidence(tmpdir):
-    grid = RasterModelGrid((3, 5))
-    grid.add_full("bedrock_surface__elevation", 0.0, at="node")
-    grid.add_full("topographic__elevation", 0.0, at="node")
+    grid = SequenceModelGrid(5, spacing=1.0)
 
     with tmpdir.as_cwd():
         with open("subsidence.csv", "w") as fp:
             fp.write("0.0,0.0\n4.0,40.0")
         subsidence = SubsidenceTimeSeries(grid, filepath="subsidence.csv")
-        assert_array_equal(
-            grid.at_node["bedrock_surface__increment_of_elevation"].reshape((3, 5)),
-            [
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-            ],
-        )
+        assert_array_equal(subsidence.subsidence_rate, [0.0, 10.0, 20.0, 30.0, 40.0])
+        assert_array_equal(grid.at_node["bedrock_surface__increment_of_elevation"], 0.0)
 
+        grid.at_node["bedrock_surface__increment_of_elevation"][:] = 0.0
         subsidence.run_one_step(dt=1.0)
         assert_array_equal(
-            grid.at_node["bedrock_surface__elevation"].reshape((3, 5)),
-            [
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-                [0.0, 10.0, 20.0, 30.0, 40.0],
-            ],
+            grid.get_profile("bedrock_surface__increment_of_elevation"),
+            [0.0, 10.0, 20.0, 30.0, 40.0],
         )
 
+
+def test_add_to_existing_subsidence(tmpdir):
+    grid = SequenceModelGrid(5, spacing=1.0)
+    grid.add_ones("bedrock_surface__increment_of_elevation", at="node")
+
+    with tmpdir.as_cwd():
+        with open("subsidence.csv", "w") as fp:
+            fp.write("0.0,0.0\n4.0,40.0")
+        subsidence = SubsidenceTimeSeries(grid, filepath="subsidence.csv")
+        assert_array_equal(subsidence.subsidence_rate, [0.0, 10.0, 20.0, 30.0, 40.0])
+        assert_array_equal(grid.at_node["bedrock_surface__increment_of_elevation"], 1.0)
+
         subsidence.run_one_step(dt=1.0)
         assert_array_equal(
-            grid.at_node["bedrock_surface__elevation"].reshape((3, 5)),
-            [
-                [0.0, 20.0, 40.0, 60.0, 80.0],
-                [0.0, 20.0, 40.0, 60.0, 80.0],
-                [0.0, 20.0, 40.0, 60.0, 80.0],
-            ],
+            grid.get_profile("bedrock_surface__increment_of_elevation"),
+            [1.0, 11.0, 21.0, 31.0, 41.0],
         )
 
 
 def test_set_subsidence_file(tmpdir):
-    grid = RasterModelGrid((3, 5))
-    grid.add_full("bedrock_surface__elevation", 0.0, at="node")
-    grid.add_full("topographic__elevation", 0.0, at="node")
+    grid = SequenceModelGrid(5, spacing=1.0)
 
     with tmpdir.as_cwd():
         with open("subsidence-0.csv", "w") as fp:
             fp.write("0.0,10.0\n4.0,10.0")
+        with open("subsidence-1.csv", "w") as fp:
+            fp.write("0.0,-100.0\n4.0,-100.0")
+
         subsidence = SubsidenceTimeSeries(grid, filepath="subsidence-0.csv")
-        assert_array_equal(
-            grid.at_node["bedrock_surface__increment_of_elevation"], 10.0
-        )
+        assert_array_equal(subsidence.subsidence_rate, 10.0)
 
         subsidence.run_one_step(dt=1.0)
-        assert_array_equal(grid.at_node["bedrock_surface__elevation"], 10.0)
-
-        with open("subsidence-1.csv", "w") as fp:
-            fp.write("0.0,-100.0\n4.0,-100.0")
-        subsidence.filepath = "subsidence-1.csv"
         assert_array_equal(
-            grid.at_node["bedrock_surface__increment_of_elevation"], -100.0
+            grid.get_profile("bedrock_surface__increment_of_elevation"), 10.0
         )
 
+        grid.at_node["bedrock_surface__increment_of_elevation"][:] = 0.0
+        subsidence.filepath = "subsidence-1.csv"
+        assert_array_equal(subsidence.subsidence_rate, -100.0)
+
         subsidence.run_one_step(dt=1.0)
-        assert_array_equal(grid.at_node["bedrock_surface__elevation"], -90.0)
+        assert_array_equal(
+            grid.get_profile("bedrock_surface__increment_of_elevation"), -100
+        )