diff --git a/tests/data/lidar_10m_tests.tif b/tests/data/lidar_10m_tests.tif
new file mode 100644
index 0000000..89e19e5
Binary files /dev/null and b/tests/data/lidar_10m_tests.tif differ
diff --git a/tests/test_intertidal.py b/tests/test_intertidal.py
index 0ef60ee..3ecc3cc 100644
--- a/tests/test_intertidal.py
+++ b/tests/test_intertidal.py
@@ -1,15 +1,25 @@
 import pytest
+import pickle
+import rioxarray
 from click.testing import CliRunner
+
 from intertidal.elevation import intertidal_cli, elevation
-import pickle
+from intertidal.validation import eval_metrics, map_raster, preprocess_validation
+from intertidal.extents import load_reproject
 
 
 @pytest.fixture()
 def satellite_ds():
+    """
+    Loads a timeseries of satellite data from a .pickle file.
+    TODO: Replace this with data loaded directly from datacube
+    after adding access to prod database.
+    """
     with open("tests/data/satellite_ds.pickle", "rb") as handle:
         return pickle.load(handle)
 
 
+@pytest.mark.dependency()
 def test_intertidal_cli():
     runner = CliRunner()
     result = runner.invoke(
@@ -28,6 +38,41 @@ def test_intertidal_cli():
     assert result.exit_code == 0
 
 
+@pytest.mark.dependency(depends=["test_intertidal_cli"])
+def test_dem_accuracy(
+    val_path="tests/data/lidar_10m_tests.tif",
+    mod_path="data/interim/testing/2020-2022/testing_2020_2022_elevation.tif",
+):
+    """
+    Compares elevation outputs of the previous CLI step against
+    validation data, and calculates and evaluates a range of accuracy 
+    metrics.
+    """
+    
+    # Load validation data
+    validation_da = rioxarray.open_rasterio(val_path, masked=True).squeeze("band")
+
+    # Load modelled elevation and uncertainty data
+    modelled_da = load_reproject(
+        path=mod_path,
+        gbox=validation_ds.odc.geobox,
+        name="modelled_ds",
+        resampling="average",
+    ).band_data
+
+    # Preprocess and calculate accuracy statistics
+    validation_z, modelled_z = preprocess_validation(modelled_ds, validation_ds, None)
+    accuracy_df = eval_metrics(x=validation_z, y=modelled_z, round=3)
+
+    # Assert accuracy is within tolerance
+    assert accuracy_df.Correlation > 0.9
+    assert accuracy_df.RMSE < 0.25
+    assert accuracy_df.MAE < 0.2
+    assert accuracy_df["R-squared"] > 0.8
+    assert accuracy_df.Bias < 0.2
+    assert abs(accuracy_df["Regression slope"] - 1) < 0.1
+
+
 def test_elevation(satellite_ds):
     ds, ds_aux, tide_m = elevation(
         satellite_ds,
@@ -44,10 +89,13 @@ def test_elevation(satellite_ds):
         study_area=None,
         log=None,
     )
+    """
+    Verify that elevation code produces expected outputs.
+    """    
 
     # Verify that ds contains correct variables
     assert "elevation" in ds.data_vars
     assert "elevation_uncertainty" in ds.data_vars
 
-    # Verify that ds is a single layer
-    assert "time" not in ds.dims
\ No newline at end of file
+    # Verify that ds is a single layer with no time dimension
+    assert "time" not in ds.dims